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1.        INTRODUCTION 

The  advent  of  space-based  weapon  systems  in  our  times  has  raised  the  prospects  of  future 
"Star  Wars"  conflicts,  rendering  the  potential  use  of  explosive  devices  against  space  targets  a  present 
day  engineering  reality.  The  warhead  of  choice  in  space  seems  to  be  of  the  fragmentation  type,  for 
obvious  reasons.  The  effectiveness  of  fragments  is  unhampered  by  the  space  environment  (lack  of  air 
may  even  be  helpful).  By  contrast,  bare  charges  in  space  are  considerably  less  efficient  than  in  air. 
One  may  wonder  why  this  is  so  since  in  air,  as  in  space,  the  same  amount  of  chemical  energy  is 
released  through  the  detonation  process.  The  explanation  is  that  the  difference  is  in  the  much  larger 
mass  involved  in  the  air  blast,  relative  to  the  bare  charge  mass. 

For  a  more  comprehensive  explanation,  we  take  a  close  look  at  the  process  by  which  an  explosive- 
driven  air  blast  wave  is  generated.  The  explosive  products  effectively  constitute  a  rapidly  expanding 
spherical  piston  (typical  initial  speed  around  6  km' sec),  which  drives  an  intense  shock  wave  into  the 
surrounding  air.  At  a  typical  range  of  100RQ  (and  with  air  density  equal  to  about  1/1000  of  charge 
density),  the  mass  of  air  entrained  by  the  shock  is  about  1000  times  the  charge  mass.  Thus,  the 
highly  concentrated  initial  explosive  energy,  has  spread  over  a  much  larger  mass  than  that  of  the 
charge,  via  the  mechanism  of  wave  propagation  in  compressible  media,  resulting  in  an  increased 
momentum.    For  a  comprehensive  treatment  of  blast  waves  in  air  the  reader  is  referred  to  Baker[l]. 

It  is  also  worthwhile  noting  that  explosive  products  in  space  typically  attain  hypersonic  speed  prior 
to  impacting  at  the  target.  The  flow  velocity  in  an  air  blast  is  typically  subsonic  or  somewhat 
supersonic.  It  is  thus  expected  that  the  actual  gasdynamic  interaction  between  the  blast  How  and  a 
stationary  target,  will  be  fundamentally  different  in  these  two  cases. 

We  contend  that  blast  effects  in  space  may  still  be  of  practical  interest  for  reasons  such  as  the 

following  : 

(i)  Notwithstanding  the  poor  efficiency  of  a  bare  charge,  its  use  should  not  be  ruled  out  altogether. 
Fragments  would  contribute  to  existing  -  and  potentially  hazardous  -  population  of  space 
debris,  underlining  the  obvious  fact  that  there  is  no  absolutely  safe  standoff  distance  from  an 
isotropic  fragmentation  warhead.   A  clean  bare  charge  may  thus  be  a  reasonable  alternative. 

(ii)  Even  a  fragmentation  warhead  has  some  residual  blast  capacity,  which  has  to  be  considered 
either  as  a  factor  in  enhancing  target  damage,  or  as  a  threat  to  be  reckoned  with  in  determining 
a  safe  standoff  distance. 


The  key  idea  of  the  present  model  is  a  combination  of  the  assumption  that  target  dynamic 
response  is  related  primarily  to  total  blast  impulse,  and  the  physically  plausible  notion  that  this 
impulse  is  equal  to  the  total  momentum  of  that  portion  of  the  expanding  explosive  products  which 
impacts  at  the  target.  The  sense  in  which  this  simple  notion  constitutes  an  approximation  to  a  proper 
gasdynamic  analysis  of  the  interaction  between  the  fluid  and  the  target,  is  clarified  in  Ch.  2.  In  that 
chapter  we  also  present  an  illuminating  comparison  between  impulsive  blast  loading  in  air  and  in 
space. 

In  order  to  demonstrate  the  ChargeMass-Range-Damage  relationship  implied  by  our  impact  blast 
approximation,  we  chose  a  simple  target  model :  A  cantilever  beam  with  a  rigid-perfectly  plastic 
stress-strain  relationship.  It  represents  an  extended  structural  element  such  as  a  solar  panel  or  an 
antenna.  We  make  use  of  studies  conducted  by  Mentel  [2]  and  by  Bodner  and  Symonds  [3],  which 
showed  that  by  and  large,  the  effect  of  accelerating  the  beam  impulsively  was  to  cause  a  rotation 
about  a  plastic  hinge  at  the  point  of  support.  The  final  angle  of  rotation  is  generally  proportional  to 
the  initial  kinetic  energy,  so  that  equating  damage  with  that  angle,  results  in  damage  being 
proportional  to  the  square  of  the  impulse  imparted  to  the  target  by  the  blast  loading.  A  presentation 
of  this  dynamic  response  model,  including  a  sample  case,  is  given  in  Ch.  3. 

Our  ChargeMass-Range-Damage  relationship  may  imply  some  far-reaching  conclusions  when 
applied  to  the  analysis  of  a  more  general  configuration  than  the  single-charge/ single-target  case.  In 
Ch.  4  we  present  a  simple  analysis  of  a  sub-munition  configuration  of  N  bare  charges,  concluding 
that  it  seems  to  have  no  advantage  in  efficiency,  relative  to  a  single  charge  of  equal  mass.  Sections  5 
and  6  contain  conclusions  and  references,  correspondingly. 

We  conclude  the  introduction  by  listing  the  main  assumptions  made  in  the  present  study  : 

(a)  Blast  loading  and  target  response  are  uncoupled.    This  is  true  since  typically  the  target  mass  is 
much  larger  than  the  mass  of  that  portion  of  the  explosive  products  which  impacts  on  it. 

(b)  Dynamic  target  response  is  independent  of  specific  loading  time  history.    It  depends  solely  on 
total  (time-integrated)  impulse. 

(c)  The  target  is  a  panel  extended  as  a  relatively  supple  cantilever.    It  is  supported  by  a  relatively 
rigid  and  massive  core  structure. 

(d)  The  charge  is  a  sphere  detonated  at  its  center.   The  expansion  is  spherically  symmetric. 

(e)  Target  surface  is  normal  to  local  flow  vector. 

(0      Target  orbital  velocity  relative  to  the  center  of  the  charge  is  negligible,  compared  with  the 
velocity  of  the  expanding  products. 


2.        IMPACT  BLAST  LOADING 

Consider  the  expanding  explosive  products  impacting  at  a  target  as  shown  in  Fig.  2-1.  By 
regarding  the  fluid  as  an  ensemble  of  non-interacting  particles  moving  at  velocity  U(R,t)  ,  and  by 
assuming  a  no-rebound  normal  impact  at  the  surface,  the  pressure  time  history  is  given  by  : 

P,(t)  =  p(R,t)[U(R,t)]2  (2-1) 


How  is  this  simple  impact  mechanism  related  to  the  actual  gasdynamic  interaction  between  the 
expanding  explosive  products  and  the  target?  When  a  target  is  located  at  a  range  of  at  least  several 
charge  radii,  two  features  in  the  free  stream  of  the  oncoming  fluid  are  significant  :  The  flow  is  highly 
hypersonic  (Mach  number  20  or  higher),  and  the  static  pressure  is  very  small,  which  means  that 
P  +  pU2  *  pU2  .  These  facts  were  born  out  by  a  numerical  computation  which  we  performed  for  a 
typical  high  explosive  characterized  by  the  following  parameters  : 

p0  =  1800    (kg  m-3) 

YCJ  "  3 

(2-2) 

DCJ  =  8    (m  ms"1) 

Q0  =  DCJ2/[2(yCJ2  -   1)]  =  4    (MJ  kg"1) 

Where  Q0  was  determined  by  assuming  that  the  detonation  corresponded  to  the  CJ  point  on  the 
explosive  Hugoniot  curve,  and  that  the  detonation  products  were  an  ideal  gas  with  a  specific-heat 
ratio  YCJ  .  The  spherically  expanding  flow  was  computed  by  integrating  the  Euler  equations  for 
isentropic  flow  via  a  high-resolution  conservative  finite-difference  scheme  [4-6].  The  initial  conditions 
were  the  self-similar  flow  field  of  a  just-detonated  spherical  charge  given  by  Taylor  [7].  The  code 
GRP  with  which  the  computation  was  performed  is  described  and  listed  in  Appendix  A. 

Consider  the  flow  at  a  stationary  target,  which  begins  at  the  moment  of  arrival  of  the  expanding 
explosive  products  (Fig.  2-2).  A  qualitative  description  of  the  ensuing  flow  pattern  is  made  by 
observing  its  evolution  in  time.  Immediately  following  the  initial  (normal)  impact,  the  fluid  is  stopped 
at  the  target  by  a  backward-propagating  shock  wave  reflected  from  the  surface.    Since  the  target  is  of 


finite  extent,  the  fluid  between  the  shock  and  the  surface  is  accelerated  laterally,  and  streamlines  that 
tend  to  curve  around  the  target  are  being  formed.  If  the  oncoming  flow  were  stationary,  the  flow 
field  would  evolve  toward  the  familiar  configuration  of  a  detached  bow-shock,  positioned  at  a 
relatively  narrow  standoff  distance  from  the  surface. 

Let  us  find  the  post-shock  pressure  in  these  two  limiting  phases.  In  the  initial  phase,  the  fluid  is 
stopped  at  the  target  by  a  reflected  shock  (Fig.  2-3a),  and  in  the  pseudo-stationary  phase  (Fig.  2-3b), 
the  shock  is  stationary.  In  either  case  we  find  the  post-shock  pressure  to  be  given  by  a  pressure- 
recovery  expression  of  the  form  : 

P2  =  apU2  (2-3) 

Where  a  is  a  constant  related  to  the  appropriate  y  (assuming  the  expanded  explosive  products 
are  an  ideal  gas).   The  governing  equations  in  the  reflected  shock  case  are  : 

p(U  +  S)  =  p2s 

p(U  +  S)2=P2  (2-4) 

P(Y  +  1)/(Y  -  1)  =  p2  (strong  shock) 

Where  the  unknowns  are  p,  ,   P,  ,  S  . 
The  equations  for  the  stationary  shock  case  are  : 

pU  =  p2u2 

pU2  =  P2  +  p2U22  (2-5) 

P(Y  +  1)/(Y  ~  1)  =  P2  (strong  shock) 

Where  the  unknowns  are  p2  ,  U2  ,  P2  .  Thus,  solving  for  a  in  the  two  cases  represented  by 
equations  (2-4)  and  (2-5),  we  get  : 


CO  00 

I(R)  =jPs(t)dt    =Jp(R,t)[U(R,t)]2dt 


Reflected  shock        a  =  [(y  +  l)/2]2 

(2-6) 
Stationary  shock      a  =  2/(y  +  l) 

In  either  case,  since  the  gas  is  not  dense,  the  effective  range  of  y  is  somewhere  between  1.0  and 
1.4  ,  so  that  setting  a  =  1  is  an  approximation  commensurate  with  the  overall  crudeness  of  the 
present  impact  blast  model.  Since  the  flow  in  the  layer  between  the  shock  and  the  target  is  low 
subsonic  (at  least  it  is  so  away  from  target  edges),  the  post-shock  pressure  is  a  reasonable  substitute 
for  the  surface  pressure.  Also,  a  =  1  is  an  appropriate  approximation  where  the  flow  is  so 
rarefied  that  it  is  collisionless.  In  this  limit,  a  =  1  corresponds  to  full  thermal  accommodation  of 
re-emitted  molecules  from  a  presumably  cold  surface. 

The  foregoing  analysis  constitutes  a  justification  of  the  impact  approximation  to  the  surface 
pressure  (2-1).  Now  we  turn  to  the  task  of  evaluating  the  impulse  which  is  defined  as  the  time- 
integrated  surface  pressure.   Using  the  impact  approximation  (2-1),  the  impulse  is  given  by  : 

to  CO 

o  0 

Let  us  introduce  a  Lagrange  mass  coordinate  m  which  enables  a  transformation  from  the  Euler 
system  (R,t)  to  the  Lagrange  system  (m,t).  The  differential  relation  associated  with  this 
transformation  at  constant   R  is  : 

dm  =  47iR2p(R,t)U(R,t)dt  (2-8) 

Since  it  is  assumed  that  the  fluid  is  not  accelerated  at  any  (R,t)  in  the  range  of  interest  for  blast 
loading,  the  velocity  U(R,t)  can  be  regarded  as  function  solely  of  the  mass  coordinate  ,  so 
that  U(R,t)  =  U(m)  .  Using  (2-8)  we  are  then  able  to  cast  the  impact  blast  expression  (2-7)  in  the 
following  simple  and  physically  appealing  form  : 

I(R)  =  Z/47TR2 

*  (2-9) 

Z  =   lu(m)dm 
o 

The  upper  limit  VV  in  (2-9),  which  is  consistent  with  the  upper  limit  °o  in  (2-7),  implies  that  the 
total  impulse  is  somewhat  overestimated,  since  it  contains  contributions  from  the  innermost  layers  of 
the  explosive  products  that  will  arrive  at  the  target  as  t  -*  °o  . 


The  total  momentum  Z  is  thus  a  constant  which  can  be  evaluated  for  any  specific  explosive 
charge  by  numerical  integration.  We  performed  this  computation  with  the  code  GRP  described  in 
Appendix  A.  In  doing  so  for  the  typical  explosive  (2-2),  we  found  out  that  the  impulse  (2-9)  was  a 
reasonable  approximation  at  ranges  as  low  as  R  =  3RQ  .  Furthermore,  it  was  found  that  Z  could 
be  approximated  by  the  maximum  attainable  momentum  for  the  given  charge  mass  and 
energy  W(2Q0)^2  ,  to  within  about  6%  .  Apparently,  the  total  momentum  is  not  overly  sensitive 
to  the  exact  velocity  distribution  function  U(m)  ,  so  that  assuming  a  value  of  Z  appropriate  to  the 
uniform  distribution  U(m)  =  (2Q0)1/2  is  a  reasonable  approximation.  Thus  we  finally  arrive  at  the 
following  closed-form  approximation  for  the  blast  impulse  : 

I(R)  =  KW(2Q0)1/2/47tR2 

(2-10) 

K  =    1 

Where  the  coefficient  K  is  retained  in  order  to  suggest  that  its  value  be  determined  more  accurately 
from  detailed  experimental  or  computational  data,  in  the  event  that  such  data  become  available.  At 
present  our  best  estimate  is  K  =  1. 

There  is  one  comparison,  however,  which  can  readily  be  made  with  available  data.  We  refer  to 
impulsive  blast  loading  in  air,  such  as  given  by  Baker  (Ref.  1,  Fig.  6.3  in  the  supplement).  The 
comparison  is  conveniently  made  with  a  non-dimensional  form  of  (2-10),  which  is  rewritten  as  : 

I  =  I(R)  [47tR02/W(2Q0)1/2]  =  (R/R0)2  (2-11) 

The  air  blast  data  has  to  be  converted  to  the  same  normalization  scheme  as  in  Eq.  (2-11),  before 
the  comparison  can  be  made.  Considering  the  definition  of  I  in  (2-11)  above,  and  the  definition  of 
scaled  range  and  air  blast  impulse  (Table  6.2  of  Ref.  I),  this  conversion  is  done  by  multiplying  the 
scaled  air  impulse  and  range  by  the  following  coefficients  (sea-level  air  is  assumed)  : 

Impulse  Multiplier  p  =  3(2y)"1/2(471/3)1/3  (Pa/PaQ0)1/6  (Pa/P0),/2  =  -01204 

Range  Multiplier  6  =  (4K/3)1!3  (P0Q0/Pa)1/3  =  67.06  (2-12) 

pa  =  1.3    (kgnf3)  Pa  =  0.1    (MPa)  y  =  1.4 


The  air  blast  conversion  was  done  by  a  small  code  which  is  given  in  Appendix  B.  The  air  and 
space  blast  impulses  are  shown  in  Fig.  2-4.  We  note  that  at  ranges  larger  than  about  10  charge  radii, 
the  air  blast  impulse  is  higher  than  the  space  impulse,  and  the  gap  widens  as  the  range  increases. 
This  observation  is  consistent  with  the  qualitative  explanation  given  in  the  introduction,  which 
attributed  this  effect  to  the  increase  in  the  entrained  air  mass  at  higher  range.  At  ranges  lower  than 
10  charge  radii,  the  air  mass  is  relatively  insignificant,  so  that  one  may  expect  the  blast  impulses  in  air 
and  in  space  to  be  comparable.  Indeed,  the  inverse-square  variation  of  impulse  with  range  is 
apparent  for  the  air  blast  at  low  range.  In  absolute  values,  however,  the  low-range  space  impulse  is 
higher  by  a  factor  of  about  1.7.  This  might  be  interpreted  as  indicating  that  choosing  k  =  1/1.7 
would  be  the  appropriate  "calibration".  However,  we  do  not  propose  to  do  so,  since  we  are  not  able 
to  trace  the  various  factors  affecting  the  low- range  impulse  as  given  by  Baker  [1];  they  may  somehow 
depend  on  the  presence  of  air,  as  well  as  on  other  parameters  such  as  target  size  and  equation  of  state 
of  the  explosion  products. 
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Figure   2-1.        Impact   Blast   Loading 
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Figure  2-2.        Shock  Reflection  at  Impact  Phase 
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Figure   2-3.        Limiting   Cases   of  Shock   Reflection 
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Figure   2-4.        Impulse   of  Normally   Reflected   Blast   Wave  at   Sea-Level  and  in   Space 


3.       TARGET  DYNAMIC   RESPONSE 

For  the  sake  of  constructing  representative  ChargeMass- Range-Damage  relations  from  our  impact 
approximation  to  the  blast  impulse  (2-10),  we  suggest  a  simple  idealized  structure  as  target  model.  It 
is  a  cantilever  beam  made  of  a  metal  characterized  by  a  rigid-perfectly  plastic  stress-strain  relation. 

This  model  is  supposed  to  represent  an  extended  spacecraft  component  such  as  a  solar  panel  or  an 
antenna.  The  core  structure  is  assumed  to  be  much  more  massive  and  rigid  than  the  extended 
structural  element,  so  that  the  cantilever  can  be  idealized  as  being  rigidly  supported.  The  sole 
dynamic  and  structural  parameters  are  hence  those  of  the  cantilever. 

For  this  purpose  we  make  use  of  an  experimental  and  theoretical  investigation  of  uniform 
cantilever  beams  subjected  to  impulsive  loading  that  was  conducted  by  Mentel  [2].  Aluminum  alloy 
beams  were  held  in  a  massive  support  that  was  gliding  along  a  rail  at  speed  V  ,  until  it  was  abruptly 
stopped  by  a  very  massive  anvil.  After  the  system  came  to  rest,  the  beams  were  observed  to  have 
rotated  through  an  angle  G   about  the  point  of  support,  with  little  deformation  elsewhere  (Fig.  3-1). 

The  theoretical  model  suggested  by  Mentel  [2]  for  predicting  0(V)  ,  can  be  described  as 
comprising  two  stages.  Immediately  following  the  impact,  the  beam  commences  rotating  rigidly 
about  the  support  point,  with  an  angular  momentum  equal  to  the  pre-collision  moment  of 
momentum  about  that  point.  This  application  of  the  principle  of  conservation  of  moment  of 
momentum  entails  an  abrupt  re-distribution  of  velocity  in  the  beam,  with  velocity  being  proportional 
to  distance  from  support,  and  the  tip  moving  at  1.5  V  .  The  angle  0  is  subsequently  determined 
from  the  requirement  that  the  rotational  kinetic  energy  be  dissipated  as  plastic  hinge  work  M  9  . 
The  resulting  0(V)   expression  is  : 

0  =  (3/8)nLV2/Mp  (3-1) 

We  now  make  one  more  step  in  formulating  the  model,  in  that  we  postulate  that  the  angle  0  is 
a  measure  of  damage  .   Using  the  following  expressions  for  M    ,  ji   and  V: 

Mp  =  (l/4)Yh2 

H  =  Pph  (3-2) 

V  =  I(R)/n 
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We  get  from  (2-10)  and  (3-1)  the  following  ChargeMass-Range-Damage  (W-R-9)  relationship  : 

r  =  CW1/2 

(3-3) 
C  =  [(3/167r2e)(LQ0/ppYh3)]1/4 

We  note  that  the  effective  range  for  a  specified  target  and  "damage  level"  0  ,   is  proportional  to 
the  square  root  of  the  charge  mass  W  . 

Using  the  data  for  the  typical  explosive  (2-2),  and  the  following  data  for  a  specific  aluminum  beam, 
we  get  for  this  sample  case  : 

h  =  0.002    (m) 

L  =  1.0    (m) 

pp  =  2700    (kg  m"3)  (3-4) 

Y  =  300    (MPa) 

C  =   1.85  6-1/4    (m  kg"1/2) 

The  ChargeMass-Range-Damage  relationship  corresponding  to  this  sample  case  is  depicted  in 
Fig.  3-2. 
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4.        CLUSTER  CONFIGURATION 

In  a  cluster  configuration,  the  gain  in  damage  is  presumably  a  result  of  a  favorable  design  tradeoff 
between  reduced  charge  mass  and  reduced  range.  Can  such  a  gain  be  achieved  for  a  space  system, 
assuming  the  ChargeMass-Range-Damage  law  (3-3)  to  hold?  It  can  be  shown  that  by  adopting 
some  simple  strategy  of  sub-munition  dispersion  and  initiation,  equation  (3-3)  implies  no  gain  in 
target  damage. 

Let  us  assume  for  the  sake  of  a  reasonably  simple  analysis,  that  dispersion  and  initiation  of  sub- 
charges  would  take  place  according  to  the  following  scheme  : 


(a)  The  N  sub-charges  appear  to  fan  out  from  a  common  virtual  center,  moving  at  equal  speeds. 
At  subsequent  times,  their  centers  are  uniformly  distributed  over  an  expanding  spherical 
envelop. 

(b)  The  target  moves  at  a  constant  velocity  relative  to  the  virtual  center.  Its  point  of  closest 
approach  to  that  center  is  at  range   R  . 

(c)  The  timing  for  dispersion  is  chosen  so  that  the  target  intersects  (tangentially)  with  the  spherical 
envelop  at  the  point  of  closest  approach  (Fig.  4-1).  This  is  also  the  point  at  which  the  blast 
from  a  single-charge  configuration  detonated  at  the  virtual  center,  would  have  impacted  at  the 
target. 

(d)  All  sub-munitions  are  detonated  at  this  "moment  of  closest  approach". 

(e)  It  is  assumed  that  each  spherical  cap  of  area  4ttR2,N  will  contain  one,  and  only  one.  sub- 
charge.  The  probability  of  the  charge  location  on  that  cap  is  assumed  to  be  uniformly 
distributed.  The  expected  location  on  the  cap  is  hence  that  latitude  line  <p  which  divides  the 
cap  into  two  parts  of  equal  area  (Fig.  4-2). 

(f)  It  is  assumed  that  the  target  is  subjected  to  the  blast  of  a  single  sub-charge,  which  is  located  on 
the  mid-area  latitude  <p    of  the  spherical  cap  that  surrounds  the  target  (Fig.  4-2). 

Since  the  area  of  the  spherical  cap  subtended  by  (p   is  4ttR"/(2N)  ,   the  angle  tp   is  given  by  : 

sin((p/2)  =  (2N)-'/2  (4-1) 

We  seek  a  comparison  between  the  deflection  6  for  a  single  charge  (W,R)  ,  and  the  deflection 
Qy  in  the  sub-munition  case  (  W\,  =  W/N  ,  RN  =  2Rsin((p/2)  ).  From  the  ChargeMass-Range- 
Damage  law  (3-3),  using  also  Eq.  (4-1),  we  get  : 
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(8N/6)  -    (WN/W)^  (R/RN)4  =  1/4 


(4-2) 


Consequently,  there  is  no  potential  gain  in  a  tradeoff  between  charge  mass  and  range,  for  a  cluster 
configuration  with  the  aforementioned  dispersion  scheme.  The  factor  1/4  ,  along  with  the  mass 
overhead  inherent  in  constructing  a  multi-charge  configuration,  indicate  that  in  causing  blast  damage, 
a  single  charge  is  more  effective  than  an  equal-mass  isotropically  dispersed  cluster. 
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Figure  4-1.        Target   Intercept  at  Closest  Approach 
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5.        DISCUSSION   AND   CONCLUSIONS 

Our  analysis  pertains  to  a  bare  explosive  charge  initiated  at  a  point  of  closest  approach  to  the 
target.  We  have  shown  that  the  loading  impulse  on  a  planform  target  is  given  by  the  impact 
approximation  (2-7),  which  states  that  the  impulse  is  proportional  to  the  charge  mass  and  inversely 
proportional  to  the  range  squared.  The  impulse  in  space  has  been  compared  with  impulse  in  air  at 
sea-level.  It  was  found  that  the  two  are  quite  comparable  at  close  range  (10  charge  radii  or  less), 
exhibiting  identical  variation  with  range.  At  far  ranges,  the  impulse  in  air  is  the  higher  one.  This  is 
consistent  with  the  notion  that  spreading  the  explosive  energy  over  larger  air  mass  results  in  larger 
momentum  (and  hence  reflected  impulse).  We  then  proceeded  to  develop  the  ChargeMass-Range- 
Damage  law  (3-3)  for  an  impulse-responsive  target,  which  states  that  blast  damage  is  proportional  to 
the  square  of  the  charge  mass  and  inversely  proportional  to  the  fourth  power  of  the  range.  These 
results  were  obtained  by  introducing  extensive  simplifications  in  the  analysis  of  gasdynamic 
interaction,  and  in  the  analysis  of  dynamic  target  response.  We  have  further  shown  that  this  damage 
law  also  implies  that  no  gain  can  be  achieved  by  an  idealized  cluster  configuration  of  bare  sub- 
charges,  relative  to  a  single  charge  of  equal  total  mass. 

It  is  worthwhile  noting  that  all  assumptions  introduced  in  the  course  of  formulating  the  impact 
blast  approximation  and  the  structural  dynamic  response  to  impulsive  loading,  imply  that  target 
damage  is  overestimated.  The  only  exception  is  the  approximation  in  setting  a  =  1  ,  which  can  be 
readily  rectified  by  assigning  to  a  the  reflected  shock,  value  given  in  (2-6).  Furthermore,  we  assumed 
that  the  pressure  at  the  midpoint  of  the  target,  is  the  pressure  everywhere  on  the  target.  Due  to  flow 
around  the  edges,  the  average  pressure  is  lower  than  the  midpoint  pressure.  Also,  targets  are  not 
everywhere  normal  to  the  flow  (and  charge/target  attitude  is  not  a  design  parameter).  Oblique  impact 
obviously  entails  reduced  target  loading.  In  the  area  of  structural  dynamic  response,  a  time- 
distributed  loading  function  generally  delivers  less  kinetic  energy  to  the  structure  than  an  impulsive 
loading  of  equal  total  impulse,  resulting  in  reduced  deformation  (damage).  Thus,  while  the  present 
model  may  be  regarded  as  an  over  estimate  when  applied  to  a  sure-fail  analysis,  it  is  particularly 
suitable  in  determining  a  sure-safe  range. 
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APPENDIX  A.        The   GRP  Code 

The  purpose  of  this  Appendix  is  to  provide  a  concise  description  of  the  GRP  code,  and  a  listing  of 
its  CHARGE  version.  It  is  intended  for  users  that  have  had  prior  experience  in  implementing 
schemes  for  solving  the  Euler  equation  of  compressible  flow.  The  theoretical  background  of  GRP 
schemes  constitutes  the  principles  on  which  the  code  is  founded.  Some  familiarity  (at  least)  with  this 
background,  as  given  in  References  4,  5  and  6  is  indispensable  to  any  implementation  of  GRP 
schemes.  Reference  4  is  recommended  as  an  introduction.  The  planar  GRP  scheme  is  fully  described 
in  Reference  5,  and  the  duct-flow  GRP  scheme  on  which  the  present  CHARGE  version  is  based  is 
given  in  Reference  6.  (In  CHARGE  version  the  flow  is  spherical  and  the  "duct''  area  is  set  to 
X(I)**2  ,  but  the  code  can  handle  any  area  variation  -  see  subroutines  CROSS  and  RATIO  below). 

In  GRP  schemes,  second-order  accuracy  is  achieved  by  considering  a  piecewise  linear  interpolation 
of  the  flow  in  each  cell  (Fig.  A-l),  from  which  second-order  accurate  fluxes  at  each  cell  interface  are 
evaluated  through  an  analysis  of  a  local  Generalized  Riemann  Problem  (GRP).  Briefly  stated,  the 
GRP  goes  one  step  further  than  the  Riemann  Problem  (RP),  in  that  it  seeks  (analytically)  the  first 
time-derivative  of  the  flow  that  evolves  as  the  "diaphragm"  is  removed  from  the  cell  interface,  at  the 
origin  of  the  centered  (X,T)  wave  paths  of  the  RP  solution.  The  major  computational  subroutines  are 
CYCEUL  where  the  integration  of  conservation  laws  is  performed,  RIEMAN  where  the  local 
Riemann  Problems  are  solved  by  Newton-Raphson  iterations,  MAGA  where  the  closed-form 
expressions  derived  from  the  GRP  analysis  [6]  are  used  to  compute  flow  time-derivatives  along  the 
contact  surface,  FLL'XE  where  all  the  previously  computed  information  is  used  to  extrapolate  the 
fluxes  to  mid-time-step  (T  +  DT/2)  which  constitutes  a  second-order  accurate  flux. 

The  plan  of  this  Appendix  is  as  follows.  Array  variables,  including  those  which  carry  conserved 
variables  (mass,  momentum  and  energy),  are  described  in  section  A.l.  This  is  followed  by 
descriptions  of  general  parameters  (A. 2),  labeled  COMMON  variables  (A. 3)  and  all  subroutines  (A. 4). 
We  conclude  by  giving  the  CHARGE  version  listing  (A. 5),  which  should  be  consulted  whenever  a 
reading  of  this  code  description  is  attempted. 

NOTE  :  The  present  CHARGE  version  was  implemented  in  a  GRP  code  version  that  had  been 
converted  to  treat  detonation  waves  as  chemically  reactive  compressible  flow.  However,  the 
detonation  scheme  is  effectively  neutralized  by  setting  QDET=0  (in  NETUNM).  All  variables 
pertaining  to  detonation,  such  as  arrays  Z(I),  DZ(I),  FIMZ(I),  ZMDOT(I)  and  labeled  COMMON 
variables  containing  Z  in  their  names,  should  be  ignored. 
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A.l        Array  Variables 

The  code  GRP  is  organized  so  that  all  major  subroutines  are  called  with  standard  list  of  array 
variables  which  represent  the  integration  scheme  (i.e.  the  conservation  laws),  local  Riemann  Problem 
solutions  and  second-order  accurate  fluxes.  Virtually  all  array  variables  are  initially  defined  in 
BEGIN  (initial  conditions),  and  are  subsequently  updated  at  each  time  step  in  CYCEUL.  The 
following  list  explains  the  meaning  of  these  variables.   Some  terms  used  in  the  list  are  defined  below. 

X(I)  grid  point  coordinate. 

U(I)  velocity  in  cell  I. 

P(I)  pressure  in  cell  I  (computed  from  equation  of  state). 

RO(I)  density   in   cell   I.     This   variable   is   time-integrated   according   to    the   law   of 

conservation  of  mass.   (Computed  in  CYCEUL). 
E(I)  total  energy  per  unit  volume  (including  kinetic  energy)  in  cell  I.    This  variable  is 

time-integrated  according  to  the  law  of  conservation  of  (total)  energy.   (Computed 

in  CYCEUL). 
DU(I)  velocity  difference  in  cell  I. 

DP(I)  pressure  difference  in  cell  I. 

DRO(I)  density  difference  in  cell  I. 

DG(I)  Lagrange  sound  velocity  difference  in  cell  I. 

DXSI(I)  the  Lagrange  coordinate  increment  defined  as  RO(I)*(X(I  +  l)-X(D),  for  cell  I. 

MIN(I)  inactive  in  present  version. 

US(I)  velocity    at    the    contact    surface    obtained    after    the    resolution    of  the    local 

discontinuity    at    X(I)    (Riemann    Problem    solution).     It    is    denoted    as    U      in 

References  4-6. 
PS(I)  pressure    at    the    contact    surface    obtained    after    the    resolution    of  the    local 

discontinuity   at   X(I)    (Riemann    Problem   solution).     It    is    denoted    as    P     in 

References  4-6. 
LTDOT(I)  time  derivative  of  US(I)  along  the  contact  surface.  (This  derivative  is  the  result  of 

the  GRP  analysis.    It  is  computed  in  MAGA.   See  Ref.  5  and  6). 
PIDOT(I)  time  derivative  of  PS(I)  along  the  contact  surface.    (This  derivative  is  the  result  of 

the  GRP  analysis.    It  is  computed  in  MAGA.   See  Ref.  5  and  6). 
FIMZ(I)  inactive  in  present  version. 

ZMDOT(I)  inactive  in  present  version. 
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TENA(I) 

FIRO(I) 
FIM(I) 
FIE(I) 
GIP(I) 

VOL(I) 

Z(I) 

DZ(I) 


momentum  per  unit  volume  RO(I)*U(I)  in  cell  I.   This  variable  is  time-integrated 

according  to  the  law  of  conservation  of  momentum.   (Computed  in  CYCEUL). 

mass  flux  at  point  X(I)  (second-order  accurate). 

momentum  flux  at  point  X(I)  (second-order  accurate). 

energy  flux  at  point  X(I)  (second-order  accurate). 

the  pressure  term  in  the  momentum  flux.    It  corresponds  to  G(U)  in  References  4 

and  6. 

volume  of  cell  I. 

inactive  in  present  version. 

inactive  in  present  version. 


Glossary  of  terms  used  in  the  array  variables  list  : 

Cell  I  -  the  cell  between  grid  points  X(I)  and  X(I  +  1).  All  cell  variables  are  averages  per  that 
interval. 

Difference  in  cell  I  -  the  difference  between  values  of  variable  at  cell  boundaries  X(I  +  1)  and  X(I). 
Those  values  are  obtained  from  "monotonized"  piecewise  linear  distribution  of  each  variable 
in  each  cell.   (Fig.  A-l). 

Second-order  accurate  flux  -  the  flux  time-derivative  at  point  X(I)  is  computed  from  the  time- 
derivatives  of  pressure  and  velocity  along  contact  surface  PIDOT(I)  and  LTDOT(I)  (in 
FLUXE).  Then  the  the  flux  is  extrapolated  to  the  centered  time  point  (T  +  DT/2),  using 
those  derivatives.  This  centered  value  is  the  second-order  flux  for  integrating  the 
conservation  laws  between  T  and  T-t-  DT. 
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A.2       Major  Parameters 

A  list  of  major  parameters  indicating  their  meaning  and  the  routine  in  which  they  are  defined,  is 
given  below.  Those  parameters  defmed  in  NETUNM  are  the  run  input.  There  is  no  reading  of  an 
input  file  in  this  version  of  GRP  code  (and  the  only  output  is  the  printed  output). 

L  number  of  grid  points  +  1   (main  program) 

LL  L  -  1     (MAIN  PROGRAM) 

T  time    (MAINO) 

DT  time  step    (MAINO) 

TMAX  maximum  time  (when  T.GE.TMAX  the  run  is  terminated)    (NETUNM) 

TMUD  time  for  which  next  printing  will  take  place    (NETUNM) 

DTMUD  printing  time  step    (NETUNM) 

NCYC  serial  number  of  time  step  (integration  cycles)    (MAINO) 

COLELA  switch    to    evaluate    cell    differences    by    Colella's    method    when    COLELA.NE.O 

(NETUNM) 

KEYMON  key  for  monotonization  scheme  (just  one  is  presently  provided  when   COLELA. EQ.O) 

(NETUNM) 

NCYCPR  frequency  of  line  printing  at  each  cycle  (time  step)    (NETUNM) 

STAB  CFL  stability  coefficient.    Must  be  smaller  than  1.    (NETUNM) 

DTBA  next  time  step  computed  from  stability  criterion    (CYCEUL) 

DTKOD  former  time  step    (MAINO) 

KDT  index  of  cell  where  DTBA  was  determined    (CYCEUL) 
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A.3       Labeled  COMMON  variables 

Labeled  COMMONS  are  used  primarily  to  transmit  data  to  and  from  routines  that  perform  the 
major  computational  steps  of  the  GRP  scheme,  i.e,  RIEMAN,  MAGA  and  FLUXE;  these  routines 
are  called  from  CYCEUL.  When  the  value  of  any  of  those  variables  is  needed  for  later  use,  whether 
for  updating  conservation  variables  (RO,  TENA,  E),  or  for  printing,  it  is  stored  in  the  appropriate 
array.  All  labeled  COMMON  variables  are  grouped  under  labels  that  indicate  their  role,  and  their 
names  are  also  mnemonic.  Generally,  suffix  L  means  Left  and  suffix  R  means  Right.  It  may  indicate 
sides  either  with  respect  to  a  cell  interface  X(I),  or  with  respect  to  the  contact  surface  which  separates 
the  Right-  and  Left-  propagating  waves  in  a  solution  to  the  local  Riemann  Problem.  We  indicate  by 
INPUT  variables  that  are  computed  prior  to  calling  the  subroutine,  and  by  OUTPUT  variables  whose 
value  was  computed  within  the  subroutine  and  constitutes  the  result  of  calling  that  subroutine. 


COMMON  /STEPO/      Parameters  related  to  the  local  Riemann  Problem.    This  is  the  first  step  in 
the  GRP  scheme. 

UL,  PL,  ROL,  CL,  GL,  SL    -    velocity,  pressure,  density,  sound  speed.  Lagrange  sound  speed  and 

entropy,  attributed  to  Left  side  of  cell  interface  at  point  X(I).    (INPUT) 
USTAR,  PSTAR     -     velocity    and    pressure    at    the    contact    surface    obtained    when    the    local 

discontinuity   is   resolved  (i.e.,   the   solution   to   the   local   Riemann   Problem).     The 

omission  of  L  or  R  suffix  indicates  that  P  and  U  are  continuous  across  the  contact 

surface.     (OUTPUT) 
RSTARL,  CSTARL,  GSTARL   -   density,  sound  speed  and  Lagrange  sound  speed  on  the  Left  side  of 

the  contact  surface.     (OUTPUT) 
WL   -    Lagrange  velocity  of  propagation  of  the  Left-moving  shock,  relative  to  the  fluid.     (OUTPUT) 
UW(6)   -  velocity  of  propagation  of  each  wave  front  (Fig.  A-3).  relative  to  the  inertial  system 

(X).     (OUTPUT) 
HELEML   -      logical  variable.     If  HELEML.EQ..TRUE.   the   Left-propagating  wave   is  a   shock. 

Otherwise  it  is  a  (centered)  rarefaction  wave.     (OUTPUT) 
NFLUX   -       integer  variable.    It  denotes  the  region  in  the  Riemann  solution  wave  structure,  which 

contains  the  point  X(I)  for  all  time.    Refer  to  Fig.  A-3  for  illustration.     (OUTPUT) 
LAMDAL,  RATEL,  TEMPL,  TEMPSL.  ZL,  ZSTARL   -   inactive. 
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COMMON /STEP  1/  Parameters  related  to  the  time-derivative  evaluation  of  the  GRP  scheme, 
performed  in  MAGA.  The  time-derivatives  of  P  and  U  along  the  contact  surface  are 
the  main  result  of  MAGA. 

DLTDT,  DPIDT  -   time-derivatives  of  velocity  and  pressure  along  contact  surface.    (OUTPUT) 
ASTARL  -      The  directional  derivative  of  U  along  the  fan  characteristic  at  the  trailing  characteristic 

of  the  Left  rarefaction  wave.    It  is  not  evaluated  when  the  Left  wave  is  a  shock.   (See 

References  4-6)     (OUTPUT) 
DGIDTL,  DRIDTL   -   time-derivatives  of  Lagrange  sound  speed  and  density  along  the  left  side  of 

the  contact  surface.    (OUTPUT) 
DSDAL  -       Lagrange  spatial  derivative  of  entropy  on  the  left  side  of  contact  surface,  prior  to 

removal  of  the  partition  at  X(I). 
SH,  RAT  -    the  cross-section  area  and  the  x-derivative  of  ln(SH).    They  are  user-defined  in  CROSS 

and  RATIO  respectively. 
DSDASL  -      entropy    derivative    used    in    the    special    "sonic"    case    (i.e,    when    NFLUX  =  2    or 

NFLUX  =  5).   See  References  5,6  for  details.    (OUTPUT) 
LAMDSL,  DZDAL,  BETACL,  DZDASL  -  inactive. 


COMMON  /GRADS/  Used  to  transmit  flow  gradients  (that  exist  in  fluid  prior  to  removal  of  the 
partition  at  X(I))  to  MAGA. 

DUDXIL,  DPDXIL,  DGDXIL,  DRDXIL.  DSDXIL  -  gradients  of  U,  P,  G.  RO,  S  (with  respect  to 
Lagrange  coordinate).  They  are  computed  in  CYCEUL  for  transmission  to  MAGA. 
(INPUT) 

DZDXIL   -       inactive. 


COMMON  /FI/  Used  to   return  values   of  updated  flux  and  cell-interface  variables   from 

FLUXE. 

FIH1,  FIH2,  FIH3   -   second-order  flux  of  mass,  momentum  flow  (just  RO*U**2)  and  energy.    They 

are  extrapolated  to  Half  the  time  step  T  +  DT/2.     (OUTPUT) 
GIH   -   the  value  of  P  at  T+  DT/2 
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UXN,  PXN,  GXN,  ROXN  -  values  of  U,  P,  G,  RO  extrapolated  to  New  time  T  +  DT,  at  cell- 
interface.  They  are  used  in  CYCEUL  to  get  tentative  (pre-monotonized)  new  cell 
differences.   (OUTPUT) 

ZXN,  FIH4,  ZMDOTL,  ZMDOTR  -  inactive. 
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A.4       Description  of  Subroutines 


MAIN   PROGRAM 

The  task  of  this  program  is  to  allocate  array  space  for  the  NMAT  arrays  required  by  the  present 
version  of  GRP  code.  The  length  of  each  array  is  L.  The  allocation  is  done  by  calling  MAINO.  This 
standard  calling  sequence  is  maintained  hereafter,  thus  facilitating  modifications. 

MAINO 

This  subroutine  functions  as  an  overall  organization  routine.  It  can  be  read  as  a  kind  of  flow- 
chart of  the  entire  computation.  First,  run  set-up  is  done  by  calling  once  to  NETUNM  (data)  and 
BEGIN  (initial  conditions).  Then  a  loop  over  time  steps  is  begun.  In  each  cycle  the  integration  by 
one  time  step  is  performed  by  calling  CYCEUL,  and  subsequently  boundary  conditions  are 
implemented  by  calling  SAFAE.  Whenever  T.EQ.TMUD,  results  are  printed  by  calling  PRINT  and 
TMUD  is  updated  by  adding  DTMUD. 

NETUNM 

Here  data  are  set  for  a  particular  run.  User  is  invited  to  modify  this  routine.  There  is  no  input 
file.  This  routine  is  called  just  once  from  MAINO.  Note  that  the  detonation  data  section  is  skipped 
when  QDET.EQ.O. 

BEGIN 

Initial  conditions  are  set-up  in  this  routine.  The  configuration  of  some  nominal  case  is  given  in 
present  version.  (In  CHARGE  version  it  is  the  detonated  spherical  charge,  using  the  Taylor  self 
similar  solution  as  initial  conditions).  User  is  called  to  modify  this  routine  so  as  to  generate  any  other 
desired  initial  configuration. 

TAYLOR 

The  purpose  of  this  routine,  along  with  ancillary  routines  INIDAT,  RUNGE  and  DERIV,  is  to 
compute  the  self-similar  Taylor  solution  [7]  of  a  detonated  spherical  charge,  and  implement  it  as 
initial  conditions  for  the  GRP  computation  of  the  ensuing  expansion.  TAYLOR  is  called  once  by 
BEGIN. 

The  core  of  the  solution  is  the  numerical  (Runge-Kutta)  integration  of  two  coupled  ordinary 
differential  equations.   The  integration  variable  is  PSI.   (The  flow  velocity  normalized  by  DCJ  is  given 
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by  U  =  EXP(-PSI)  ).  The  two  dependent  variables  are  X  -  the  normalized  radial  coordinate  (X=  1  at 
the  sphere  boundary),  and  C  -  the  normalized  speed  of  sound.  The  integration  is  carried  out  by 
calling  RUNGE,  which  in  turn  calls  DERIV  for  the  evaluation  of  derivatives.  Data  for  the  TAYLOR 
computation  is  set  up  by  calling  (just  once)  INIDAT. 

The  initial  conditions  needed  in  BEGIN  are  values  of  mass,  momentum  and  (total)  energy  per  cell. 
These  are  most  accurately  computed  by  spatially  integrating  the  Taylor  solution,  resulting  in  lumped 
mass,  momentum  and  energy  per  cell,  which  are  then  divided  by  the  cell  volume.  This  refinement  is 
significant  since  gradients  are  high  near  the  charge  boundary  (X=  1).  A  total  mass  and  energy  check 
for  the  entire  sphere  is  performed  and  printed. 

INIDAT,     RUNGE,    DERIV 

Subroutines  used  only  in  conjunction  with  the  Taylor  initial  conditions  setup.  See  TAYLOR 
above. 

RATIO,    CROSS 

User-defined  routines.  If  A(X)  is  the  duct  cross-section  area,  then  CROSS(X)  =  A(X)  and 
RATIO(X)=D[ln(A(X))]/DX  . 

CYCEUL 

This  is  the  central  computation  routine.  All  major  stages  of  the  GRP  scheme  are  performed  by 
calling  specific  subroutines  from  CYCEUL.  Then  RO(I),  TENA(I)  and  E(I)  are  updated  to  new  time 
T  +  DT  by  solving  the  appropriate  conservation  laws  in  CYCEUL. 

The  first  loop  (DO  1)  performs  a  set  of  preparatory'  steps  as  follows  : 

(a)  CALL  RIEMAN   -   Solving  the  local  Riemann  Problem  at  each  X(I). 

(b)  CALL  MAG  A   -    Solving  the  local  Generalized  Riemann  Problem  at  each  X(I). 

(c)  CALL  FLUXE   -   Computing  second-order  fluxes  at  X(I). 

(d)  Evaluation  of  cell-interface  finite  differences  DU(I),  DP(I),  DRO(I)  in  each  cell.    These  will  be 
used  at  the  future  time  step  (after  monotonization)  for  piecewise-linear  interpolation  of  the  flow 

in  each  ceil.    (See  definition  of  DUDXIL,  DPDXIL just  preceding  the  call  to  MAGA  in  this 

loop). 

Note  that  in  present  CHARGE  version  additional  computation  of  PRESS,  PULSE1,...,  PULSE4 
has  been  added.    It  is  just  informative  and  does  not  interfere  in  any  way  with  the  execution  of  the 
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GRP  scheme.  The  purpose  of  this  computation  is  to  monitor  the  numerical  solution  and  to  observe 
the  accuracy  within  which  the  asymptotic  value  of  the  momentum  integral  Z  (Eq.  2-9  above)  is 
approached. 

In  the  second  loop  (DO  2),  the  integration  of  the  three  conservation  laws  is  performed,  using 
second-order  fluxes  that  had  been  computed  in  loop  1.  Flow  variables  such  as  P(I)  and  U(I)  are 
computed  in  this  loop  from  the  conserved  variables.  The  cycle  computation  is  concluded  by  calling 
BDOK1  for  monotonization  of  DU(I),  DP(I)  and  DRO(I). 

SAFAE 

In  this  routine  user-defined  boundary  conditions  are  implemented.  Present  version  (CHARGE) 
contains  rigid  wall  at  the  center  of  the  sphere  X(2)  =  0,  and  an  "open  boundary"  at  the  outer 
computational  zone  limit  X(L).  The  rigid  wall  condition  is  achieved  by  setting  up  a  virtual 
antisymmetric  cell  next  to  the  boundary  cell,  so  that  the  solution  to  the  local  Riemann  Problem  will 
result  in  a  non-moving  contact  surface  (USTAR  =  0).  The  open  boundary  is  an  approximation  to  an 
ideally  non-reflecting  boundary.  Here  the  virtual  cell  is  I  =  L,  and  the  flow  in  it  is  defined  as  a 
"continuation"  of  the  flow  in  the  adjacent  last  cell  1  =  LL. 

BDOK1 

Here  the  tentative  cell-interface  differences  DV(I)  are  monotonized  according  to  neighboring 
average  cell  values  V(I-l),  V(I)  and  V(I+1).  The  basic  idea  is  that  the  cell-interface  slope  DV(I) 
should  have  the  same  sign  as  the  average  slope  V(I  +  l)-V(I-l).  When  V(I)  is  a  local  extremum  DV(I) 
is  set  to  zero.  Also,  the  absolute  value  of  DV(I)  is  constrained  so  that  the  jump  from  a  cell-interface 
value  to  the  adjacent  average  value  V(I),  will  never  be  of  opposite  sign  to  DV(I). 

DCOLE 

When  COLELA  option  is  used  (not  in  present  CHARGE  version),  the  pre-monotonized  slopes 
are  simply  the  centered  difference  (V(I+  l)-V(I-l))/2.  Note  that  even  under  this  option,  the 
monotonization  routine  BDOK1  is  subsequently  called. 

PRINT 

Printing  of  results.  Reading  this  routine  is  self-explanatory.  Note  some  features  added  for 
present  CHARGE  version.    User  is  called  to  modify  this  routine  to  his  specific  needs. 
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SOF 

Run  termination  when  an  error  has  been  detected.  I  STOP  is  an  informative  index.  All  printing 
of  relevant  information  should  be  done  at  the  calling  routine  prior  to  calling  SOF.  Note  that  the  run 
is  ended  in  SOF  by  deliberately  causing  a  system  error  of  computing  SQRT(-l).  This  is  done  in  order 
to  trigger  printing  of  the  sequence  of  calling  routines  by  the  operating  system. 


RIEMAN 

Here  a  single  Riemann  Problem  (RP)  is  solved  by  calling  RIEMAN  from  CYCEUL.  Referring 
to  Fig.  A-2,  the  RP  is  solved  by  finding  the  point  of  intersection  (USTAR,PSTAR)  of  Left- 
propagating  and  Right-propagating  shock/rarefaction  adiabats  in  the  (U,P)  plane.  Prior  to  the  actual 
computation,  the  qualitative  wave  structure  is  determined.  It  is  characterized  by  the  index  NCASE  as 
follows  : 

NCASE  =  1   -   Left  wave  is  rarefaction,  Right  wave  is  shock. 
NCASE  =  2  -   Both  waves  are  shock. 

NCASE  =  3   -   Left  wave  is  shock.  Right  wave  is  rarefaction. 
NCASE  =  4  -   Both  waves  are  rarefaction. 

The  computation  of  (USTAR.PSTAR)  is  coded  separately  for  each  case.  Newton-Raphson 
iteration  is  employed,  the  first  guess  being  the  intersection  of  the  Left  and  Right  rarefaction  branches 
(or  their  extrapolations),  which  is  done  in  closed-form.  Since  in  a  smooth  flow  this  guess  is  close  to 
the  exact  (L'STAR,PSTAR),  little  extra  CPU  effort  is  spent  on  subsequent  Newton-Raphson 
iterations.   These  are  truly  needed  only  in  regions  of  shock  wave  computation. 

The  computation  in  RIEMAN  is  concluded  by  computing  UW(1),...,UW(5)  (UW( 6)  =  infinity). 
From  these  wave  speeds,  the  flux  index  NFLUX  that  denotes  the  location  of  the  X-axis  on  the  (X,T) 
wave  diagram  of  the  RP  solution  (Fig.  A-3),  is  evaluated.    It  is  later  needed  in  subroutine  FLUXE. 
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MAGA 

The  major  purpose  of  this  routine  is  to  compute  DUIDT  and  DPIDT  along  the  contact  surface 
of  the  RP  solution.  Since  U  and  P  are  continuous  across  the  contact,  so  are  their  time-derivatives 
along  the  contact.  Thus,  DUIDT  and  DPIDT  are  solved  from  a  set  of  two  linear  equations.  The 
coefficients  of  each  equation  are  determined  by  GRP  analysis  of  the  wave  on  one  side.  See 
References  4-6  (particularly  Ref.  6)  for  details. 

FLUXE 

The  major  task  of  this  routine  is  to  compute  second-order  fluxes.  This  is  done  in  two  phases. 
The  first  phase  is  up  to  statement  9  CONTINUE,  where  using  NFLUX  the  X-suffixed  values  of  flow 
variables  and  their  time-derivatives  are  defined.  An  X-suflix  means  that  the  variable  or  its  time- 
derivative  are  related  to  the  line  X  =  X(I)  on  the  (X,T)  wave  diagram  (Fig.  A-3).  In  the  second  phase, 
these  variables  and  their  time-derivatives  are  used  to  extend  fluxes  at  X(I)  to  Half-time-step  (hence 
the  suffix  H),  i.e.  T  +  DT/2.  It  is  these  fluxes  which  are  the  second-order  accurate  fluxes  for  the 
integration  of  the  conservation  laws  from  T  to  T  +  DT.  Also,  cell-interface  flow  variables  (suffix  N) 
are  extended  to  New  time  level  T  +  DT.  These  are  later  used  in  defining  cell  differences  DU(I),  DP(I) 
and  DRO(I)  in  CYCEUL. 
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A.5       Listing  of  GRP  Code 

CSOPTIONS  LIST                                         ^ga  vettHo*            CHA0001 

IMPLICIT  REAL*8(A-H,0-Z,$)                       °"' £  "fc*s'off-    CHA0002 

C   PROGRAM  GRP  -  GENERALIZED  RIEMANN  PROBLEM.  CHA0003 

C   EXPANSION  OF  A  DETONATED  SPHERICAL  CHARGE  IN  VACUUM.  CHA0004 

C   INITIAL  CONDITIONS  FROM  TAYLOR'S  SELF  SIMILAR  SOLUTION.  CHA0005 

COMMON  B(102,26)  CHA0006 

1       ,ENDB  CHA0007 

COMMON  /AB/A(50)  CHA0008 

EQUIVALENCE  (L  ,  A(  1)  )  ,  (LL  ,  A(2)  )  ,  (T,  A(3)  )  ,  (DT,  A(4)  )  ,  (TMAX,  A(5)  )  ,     CHA00  09 

1  (TMUD,A(6)),(DTMUD,A(7)),(J0B,A(8)),(NERI,A(9)),  CHA0010 

2  (JJJ,A(10)),(KEYMON,A(11)),(NCYC,AU2))  CHAOOll 
EQUIVALENCE  (COLELA, A( 13) )  CHA0012 
EQUIVALENCE  ( LAGEUL , A( 14) )  CHA0013 
EQUIVALENCE  (UGAL,A(15))  CHA0014 
EQUIVALENCE  ( KEYEK, A( 16 ) )  CHA0015 
EQUIVALENCE  (NCYCPR, A( 17 ) )  CHA0016 
EQUIVALENCE  ( STAB, A( 18) ) , ( DTBA, A( 19) ) , (DTKOD, A(20) ) , (KDT, AC 21 ) )  CHA0017 
COMMON  /M0NIT/CASEAV(4),NC14(4),NF16(6),  CHA0018 

1                 NM0NU(4),NM0NP(4),NM0NR0(4),NM0NZ(4)  CHA0019 

DIMENSION  NZER0C26)  CHA0020 

EQUIVALENCE  ( NZERO( 1 ) , NC14( 1 ) )  CHA0021 
COMMON/PULS/PRESS(10),PULSE1(10),PULSE2(10),PULSE3(10),PULSE4(10)  CHA0022 

CXXXX»XXXXXXXXXXXXX5(XXXXXXX5(XXXXXXXXXXXXXXXXXXXXXXXXXXX5«XXXXXXXXXXXXXXXXCHA0023 

DO  20  N=l,26  CHA0024 

20  NZERO(N)=0  CHA0025 
DO  21  N  =  l,4  CHA0026 

21  CASEAV(N)=0.  CHA0027 
DO  31  N=l,10  CHA0028 
PRESS(N)=0.  CHA0029 
PULSE1(N)=0.  CHA0030 
PULSE2(N)=0.  CHA0031 
PULSE3(N)=0.  CHA0032 
PULSE4(N)=0.  CHA0033 

31    CONTINUE  CHA0034 

NMAT=26  CHA0035 

C      L=(LOCF(ENDB)-tOCF(B(l,l)))/NMAT  CHA0036 

L=102  CHA0037 

LL=L-1  CHA0038 

NN=NMAT*L  CHA0039 

DO  1  1=1, L  CHA0040 

DO  1  11=1, NMAT  CHA0041 

1     B(I,II)=0.  CHA0042 

CALL  MAIN0(L,B(1,  1),B(1,  2),B(1,  3),B(1,  4), 8(1,  5),  CHA0043 

1  B(l,  6),B(1,  7),B(1,  8),B(1,  9),B(1,10),  CHA0044 

2  B(1,11),B(1,12),B(1,13),B(1,14),B(1,15),  CHA0045 

3  B(1,16),B(1,17),B(1,18),B(1,19),B(1,20),  CHA0046 

4  B(1,21),B(1,22),B(1,23),B(1,24),B(1,25),  CHA0047 

5  B(l,26))  CHA0048 
STOP  CHA0049 

END .    CHAQ050 

SUBROUTINE   MAINO  1*Pr\H0          CHA0051 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0052 

2  US,PS,UIDOT,PIDOT,  CHA0053 
*              FIMZ,ZMDOT,  CHA0054 

3  TENA,FIRO,FIM,FIE,UlP,VOL,Z,DZ)  CHA0055 
IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0056 
DIMENSION  X(L),U(L),P(L),R0(L),G(L),E(L),DU(L),DP(L),DR0(L),  CHA0057 

1  DG(L),DXSI(L),MIN(L),  CHA0058 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  CHA0059 

3  ,TENA(L),FIR0(L),FIM(L),FIE(L)  CHA0060 

4  ,GIP(L),V0L(L),Z(L),DZ(L)  CHA0061 

5  ,FIMZ(L),ZMDOT(L)  CHA0062 
COMMON  /AB/AC50)  CHA0063 
EQUIVALENCE  ( LL , A( 2) ) , (T, A( 3 ) ) , ( DT , A( 4 ) ) , (TMAX, A( 5) ) ,  CHA0064 

1  (TMUD,A(6)),(DTMUD,A(7)),(J0B,A(8)),(NERI,A(9)),  CHA0065 

2  (JJJ,A(10)),(KEYMON,A(11)),(NCYC,A(12))  CHA0066 
EQUIVALENCE  ( LAGEUL , A( 14) )  CHA0067 
EQUIVALENCE  ( NCYCPR, A( 17 ) )  CHA0068 
EQUIVALENCE  ( STAB, A( 18 )),( DTBA, A( 19 )),( DTKOD, A(20) ),( KDT, A( 21 ) )  CHA0069 
COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  CHA0070 

CXXXXX^^XXXXXXX^X3<X5(XXXXXXX^XX^XXXX^XXXXX^XXXX3(XXXXXX)(X)(XXXXXXXXXXXXXXXXCHA0O71 

T=0.  CHA0072 
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NCYC=0  CHA0073 

JJJ=0  CHA0074 

CALL  NETUNM  CHA0075 

DELT=DT  CHA0076 

CALL  BEGIN  CHA0077 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0078 

2  US,PS,UIDOT,PIDOT,  CHA0079 
X              FIMZ,ZMDOT,  CHA0080 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0081 
CALL  SAFAE  CHA0082 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0083 

2  US,PS,UIDOT,PIDOT,  CHA0084 
X              FIMZ,ZMDOT,  CHA0085 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0086 

1  NCYC=NCYC+1  CHA0087 
C  TIME  STEP  CONTROL.  CHA0088 

DT=DTBA  CHA0089 

IFCDT.GT.1.1DO*DTKOD.AND.DTKOD.NE.O.)  DT=1 . IDOXDTKOD  CHA0090 

IFCNCYC.EQ.2)  DT=DT/10.D0  CHA0091 

IF  (NCYC.EQ.l)  DT=0.  CHA0092 

IFCDT.EQ.O.)  GO  TO  11  CHA0093 

NHAD=((TMUD-T)/DT-1.D-10)  CHA0094 

IF(NHAD.GE.IO)  GO  TO  11  CHA0095 

DT=(TMUD-T)/DFLOAT(NHAD+l)  CHA0096 

11    CONTINUE  CHA0097 

T=T+DT  CHA0098 

IFCCNCYC/NCYCPR)*NCYCPR.NE.NCYC.AND.NCYC.GT.NCYCPR)  GO  TO  33  CHA0099 

PRINT  10,  NCYC,T,DT,KDT  CHA0100 
10    F0RMATC1X,'NCYC=M4,3X,  'T=  » ,  Dll  .4,  3X,  »DT=  ' ,  Dll  .  4,  3X,  »KDT  =  M4)    CHA0101 

33    CONTINUE  CHA0102 

DTBA=DTMUD  CHA0103 

KDT=0  CHA0104 

NERI=1  CHA0105 

IF  (DABS(T-TMUD) .LT.l.D-8)  NERI=0  CHA0106 

CALL  CYCEUL  CHA0107 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0108 

2  US,PS,UIDOT,PIDOT,  CHA0109 
*                                  FIMZ,ZMDOT,  CHA0110 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0111 
CALL  SAFAE  CHA0112 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0113 

2  US,PS,UIDOT,PIDOT,  CHA0114 
X             FIMZ,ZMDOT,  CHA0115 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0116 
IF  (NERI.NE.O)  GO  TO  2  CHA0117 
CALL  PRINT  CHA0118 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0119 

2  US,PS,UIDOT,PIDOT,  CHA0120 
X             FIMZ,ZMDOT,  CHA0121 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0122 
IF  (DABS(T-TMUD). LT.l.D-8)  TMUD=TMUD+DTMUD  CHA0130 

2  CONTINUE  CHA0131 
DTKOD=DT  CHA0132 
IF  (T.LT.TMAX-l.D-8)  GO  TO  1  CHA0133 
RETURN  CHA0134 

END  , CHAQJ.55 

SUBROUTINE  NETUNM  H£Tl)N+A  CHA0136 

IMPLICIT  REALX8CA-H,0-Z,$)  CHA0137 

COMMON  /AB/AC50)  CHA0138 

EQUIVALENCE  CL,AC1))  CHA0139 

EQUIVALENCE  C LL , AC  2) ) , CT, AC  3 ) ) , C DT, AC  4) ) , CTMAX, AC5) ) ,  CHA0140 

1  CTMUD,AC6)),CDTMUD,AC7)),CJ0B,AC8)), CNERI,AC9)),  CHA0141 

2  CJJJ,AC10)),CKEYMON,AC11)),  CNCYCAC12)  )  CHA0142 
EQUIVALENCE  C COLELA, AC  13) )  CHA0143 
EQUIVALENCE  ( LAGEUL , AC 14 ) )  CHA0144 
EQUIVALENCE  C KEYEK, AC  16 ) )  CHA0145 
EQUIVALENCE  ( NCYCPR, AC  17 ) )  CHA0146 
EQUIVALENCE  C STAB, AC  18 ) ) , C DTBA, AC  19 ) ) , C DTKOD, AC  20 ) ) , CKDT, AC  21 ) )  CHA0147 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET,  CHA0148 

1  RATE,TEMPC  CHA0149 

C0MM0N/DIFFUS/U2,P2,R02,ARN  CHA0150 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX  CHA0151 
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1             ,XMIN,XMAX,SMIN,SMAX,IVERSA  CHA0152 

COMMON  /GAM/GAMA , NG, MU2 , Gl , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , Gl 0 , Gl 1  CHA 0 153 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA0154 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0155 
REALX8  NG,MU2  CHA0156 
NAMELIST  /IN/LIN,GAMA,DT,TMUD,DTMUD,TMAX,  CHA0157 

1  GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX,  CHA0158 

2  SMIN,SMAX,IVERSA,KEYMON,COLELA,STAB  CHA0159 

3  ,LAGEUL,KEYEK  CHA0160 
<+            ,QDET  CHA0161 

exxxxxxxxxxxxxxxxmxxxx3«xxxxxxxxxxxxxxxxxxxxxxxx3«5€XXXX3€5«xxxxxxxxxxxxxxxcHA0162 

LIN=L  CHA0163 

LAGEUL=2  CHA0164 

NCYCPR=1  CHA0165 

KEYEK=1  CHA0166 

TMUD=0.  CHA0167 

DTMUD=10.D0  CHA0168 

TMAX=100.D0  CHA0169 

STAB=0.5D0  CHA0170 

DT=l.D-2  CHA0171 

KEYMON=l  CHA0172 

GAMA=3.D0+l.D-6  CHA0173 

QDET  =  0.04D0  CHA0174 

RATE=0.  CHA0175 

TEMPC=1.D50  CHA0176 

GODELX=16D0  CHA0177 

GODELY=20.D0  CHA0178 

IVERSA=100  CHA0179 

UMIN=0.  CHA0180 

UMAX=  1.D0  CHA0181 

PMIN=0.  CHA0182 

PMAX=0.5D0  CHA0183 

ROMIN=0.  CHA0184 

ROMAX=3.D0  CHA0185 

SMIN=0.  CHA0186 

SMAX=0.03D0  CHA0187 

COLELA=0.  CHA0188 

READ  IN  CHA0189 

CHA0190 

PRINT  IN  CHA0191 

CHA0192 

GG=2.D0*GAMA/(GAMA-1.D0)  CHA0193 

NG=GG  CHA0194 

)   CONTINUE  CHA0195 

MU2=(GAMA-1.D0)/(GAMA+1.D0)  CHA0196 

G1=GAMA-1.D0  CHA0197 

G2=1.D0-MU2  CHA0198 

G3=2.D0/(3.D0*GAMA-1.D0)  CHA0199 

G4=(GAMA+1 .DO)/2.DO  CHA0200 

G5=0.5D0*(3.D0*GAMA-1.D0)/(GAMA+1.D0)  CHA0201 

G6=(GAMA+1 . DO)/(2.DO*GAMA)  CHA 02 02 

G7=2.D0/(GAMA-1.D0)  CHA0203 

G3=(GAMA-1.D0)/(2.D0*GAMA)  CHA0204 

G9  =  (  GAMA  +  1.  DO  )/(<*.  DO  XGAMA)  CHA  020  5 

G10=1.D0/GAMA  CHA0206 

G11=(GAMA+1.D0)/4.D0  CHA0207 

G12=GAMA/(GAMA-1.D0)  CHA0208 

G13=0.5D0*(GAMA-3.D0)/(GAMA+1  .DO)  CHA0209 

G14=0.5D0*(3.D0*GAMA-5.D0)/(GAMA+1.D0)  CHA0210 

G15=GAMA*(3.D0*GAMA-1.D0)  CHA0211 

G16=(GAMA+1 . D0)/(2.D0*(GAMA-1.D0))  CHA0212 

G17=GAMA+1 .DO  CHA0213 

G18=GAMA*(GAMA+1 . DO )/ ( 3 . D0*GAMA-1 -DO)  CHA 02 14 

G19=(3.D0*GAMA-1.D0)/(GAMA+1 .DO)  CHA 02 15 

G20=2.D0*(GAMA-1 . DO )/ ( 3 . D0*GAMA-1 . D0)**2  CHA 02 16 

G21 =GAMA*( 3 . D0*GAMA-5 . DO )/ ( 3 . DOXGAMA-1 . DO ) **2  CHA 02 17 

GODELX=GODELX/2.54D0  CHA0213 

GODELY=GODELY/2.54D0  CHA0219 

CALL  NAMPLT(IVERSA)  CHA0220 

CALL  LIMIT(IOOO.DO)  CHA0221 

CALL  PLOT(0.,0.5D0,-3)  CHA0222 

PODET=0.  CHA0223 
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RO0DET=0.  CHA0224 

PCJDET=0.  CHA0225 

UCJDET=0.  CHA0226 

DCJDET=0.  CHA0227 

RCJDET=0.  CHA0228 

IFCQDET.LE.O.)  GO  TO  100  CHA0229 

C   DETONATION  DATA  CHA0230 

QDET=0.04D0  CHA0231 

P0DET=0.  CHA0232 

RO0DET=1.8D0  CHA0233 

PCJDET=P0DET-(GAMA-1 . D0)*(-QDET)*RO0DET+  CHA0234 

1  DSQRT(((GAMA-l.D0)*QDET*RO0DET)**2-2.D0*MU2*GAMA*  CHA0235 

2  C-QDET)*P0DET*RO0DET)  CHA0236 
RCJDET=RO0DET*((GAMA+l.D0)*PCJDET-P0DET)/(GAMA*PCJDET)  CHA0237 
CCJ=DSQRT(GAMA*PCJDET/RCJDET)  CHA0238 
DCJDET=CCJ*RCJDET/RO0DET  CHA0239 
UCJDET=DCJDET-CCJ  CHA0240 
PRINT  101  CHA0241 

101  FORMATdHl,/,  IX,  'DETONATION  DATA'/)  CHA0242 
PRINT  102,  QDET,GAMA,TEMPC,RATE  CHA0243 

102  F0RMATC/1X, 'QDET, GAMA, TEMP, RATE= ' , 4D18 .8)  CHA0244 
PRINT  103,  RO0DET,P0DET  CHA0245 

103  FORMATC/1X, 'UNBURNED  STATE       ROODET, P0DET= ' , 2D18 . 8)  CHA0246 
PRINT  104,  DCJDET,PCJDET,RCJDET,UCJDET  CHA0247 

104  FORMATC/1X, 'CJ  POINT      DCJDET, PCJDET, RCJDET, UCJDET= ' , 4D18 .8)     CHA0248 
100   CONTINUE  CHA0249 

RETURN  CHA0250 

END CHA0251 

SUBROUTINE    BEGIN  RP/Tf/V       CHA0252 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  °C(J  '    CHA0253 

2  US,PS,UIDOT,PIDOT,  CHA0254 
*             FIMZ,ZMDOT,  CHA0255 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0256 
IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0257 
DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  CHA0258 

1  DG(L),DXSI(L),MIN(L),  CHA0259 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  CHA0260 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  CHA0261 

4  ,GIP(L),VOL(L),Z(L),DZ(L)  CHA0262 

5  ,FIMZ(L),ZMDOT(L)  CHA0263 
COMMON  /AB/AC50)  CHA0264 
COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  CHA026  5 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA0266 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0267 
REAL*8  NG,MU2  CHA0268 
EQUIVALENCE  (LL,A(2))  CHA0269 
EQUIVALENCE  ( LAGEUL , A( 14) )  CHA0270 
EQUIVALENCE  ( UGAL , A( 15 ) )  CHA0271 
EQUIVALENCE  (STAB, AC  18 ) ) , ( DTBA, A( 19 ) ) , ( DTKOD, A( 20 ) ) , ( KDT, A( 21) )  CHA027  2 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,PODET,RO0DET,  CHA027  3 

1             RATE,TEMPC  CHA0274 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX  CHA027  5 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA  CHA0276 
COMMON/GIT/ROLIM,ELIM,XGIT(200),ROGIT(200),ROUGIT(200),EGIT(20  0)   CHA0277 

COMMON  /GITN/NPO  CHA0278 

LOGICAL  CSOF  CHA0279 
CXXXXX^X^^XXXX^XX^XXX^XXXX^^MXX^5(^5(^^5«X^X^3(XXXXXX^^X5(XXXXXXXXX5eXX5(XX^^^^CHA0  28  0 

DTBA=0.  CHA0281 

DTKOD=0.  CHA0282 

KDT=0  CHA0283 

P0=l.D-9  CHA0284 

RH00=l.D-7  CHA0285 

UO=UCJDET  CHA0286 

UGAL=0.  CHA0287 

X0=0.  CHA0288 

X1=50.D0  CHA0289 

XCHARG=10.D0  CHA0290 

XMIN=X0  CHA0291 

XMAX=X1  CHA0292 

DX=(X1-X0)/(L-2.D0)  CHA0293 

DO  1  1=2, L  CHA0294 

X(I)=X0+(I-2.D0)*DX  CHA0295 
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CONTINUE  CHA0296 

X(L)=X1  CHA0297 

U0=U0XCROSS(XCHARG)  CHA0298 

DO  2  1=2, LL  CHA0299 

IFCI.GT.2)  U(I)=U0/CROSS(X(I))  CHA0300 

P(I)=P0  CHA0301 

RO(I)=RHOO  CHA0302 

Z(I)=0.  CHA0303 

GO  TO  (31,32),  LAGEUL  CHA0304 

11  CONTINUE  CHA0305 
E(I)=P(I)/((GAMA-l.DO)xRO(I))+0.5DOXU(I)x*2+Z(I)XQDET  CHA0306 
GO  TO  30  CHA0307 

12  CONTINUE  CHA0308 

E(I)=P(I)/(GAMA-1.DO)+0.5DO*RO(I)XU(I)**2+Z(I)*RO(I)XQDET  CHA0309 

10    CONTINUE  CHA0310 

G(I)=DSQRT(GAMA*P(I)XRO(I))  CHA0311 

1    CONTINUE  CHA0312 

DO  3  1=2, LL  CHA0313 

TENA(I)=R0(I)XU(I)  CHA0314 

V0L(I)=(X(I+1)-X(I))X(X(I+1)XX2+X(I+1)XX(I)+X(I)XX2)/3.D0  CHA0315 

I    CONTINUE  CHA0316 

CHA0317 

INSERT  DETONATED  CHARGE  FLOW  FIELD  FROM  TAYLOR'S  SOLUTION.  CHA0318 

CHA0319 

CALL  TAYLOR(GAMA)  CHA0320 

RONORM=RCJDET  CHA0321 

RUNORM=RCJDET*UCJDET  CHA0322 

EN0RM=RCJDET*DCJDET**2  CHA0323 

XLIM=XGIT(NPO)  CHA0324 

NGIT=NP0-1  CHA0325 

XG1=XLIM  CHA0326 

XG2=XGIT(NGIT)  CHA0327 

AROIP  =ROGIT  (NPO)+ROLIMXXLIMXX3/3.DO  CHA0328 

AROUIP=ROUGIT(NPO)  CHA0329 

AEIP   =EGIT   (NPO)+  ELIMXXLIMXX3/3.D0  CHA0330 

XP=X(2)/XCHARG  CHA0331 

DO  100  1=2, LL  CHA0332 

IP=I+1  CHA0333 

XI=XP  CHA0334 

AROI  =AROIP  CHA0335 

AROUI=AROUIP  CHA0336 

AEI   =AEIP  CHA0337 

XP=X(IP)/XCHARG  CHA0338 

IFCDABS(XP-l.DO) .LT.l.D-10)  XP=1.D0  CHA0339 

CS0F=(XP.GE.1.D0)  CHA0340 

IF(XP.GE.XLIM)  GO  TO  101  CHA0341 

UNIFORM  FLOW  REGION  CHA0342 

DELVOL=(XLIM-XP)X(XLIM**2+XLIM*XP+XPXX2)/3.DO  CHA0  343 

AROIP  =ROGIT  (NP0)+ROLIM*DELVOL  CHA0344 

AROUIP=ROUGIT(NPO)  CHA0345 

AEIP   =EGIT   (NPO)+  ELIMXDELVOL  CHA0346 

GO  TO  102  CHA0347 

.01   CONTINUE  CHA0348 

NON  UNIFORM  FLOW  REGION.  CHA0349 

IF( .NOT.CSOF)  GO  TO  104  CHA0350 

LAST  POINT.   (THIS  IS  THE  DETONATION  FRONT  POINT   X=l) .  CHA0351 

AROIP=  0.  CHA0352 

AROUIP=0.  CHA0353 

AEIP=   0.  CHA0354 

GO  TO  102  CHA0355 

.04   CONTINUE  CHA0356 

IF(XP.LE.XG2)  GO  TO  103  CHA0357 

NGIT=NGIT-1  CHA0358 

IF(NGIT.LE.O)  CALL  SOFCBEGIN  104.  NGIT.LE.O.')  CHA0359 

XG1=XG2  CHA0360 

XG2=XGIT(NGIT)  CHA0361 

GO  TO  104  CHA0362 

.03   CONTINUE  CHA0363 

FRAC=(XP-XG1)/(XG2-XG1)  CHA0364 

IF(FRAC.LT.O.   )  CALL  SOFCBEGIN  103.  FRAC.LT.O.')               CHA0365 

IF(FRAC.GT.l.DO)  CALL  SOFCBEGIN  103.  FRAC.GT.l.1)               CHA0366 

AROIP  =(l.D0-FRAC)*ROGIT  (NGIT+1 )+FRACXROGIT  (NGIT)  CHA0367 
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AROUIP=(l .D0-FRAC)XROUGIT(NGIT+1)+FRAC*ROUGIT(NGIT)  CHA0368 

AEIP   =(1.D0-FRAC)*EGIT   (NGIT+D  +  FRACXEGIT   (NGIT)  CHA0369 

102   CONTINUE  CHA0370 

C   COMPUTE  MASS,  MOMENTUM  AND  ENERGY  DENSITIES.  CHA0371 

IF(XP.LE.XLIM)  GO  TO  105  CHA0372 

C   CONSERVATION-FORM  DEFINITION  OF  MASS,  MOMENTUM  AND  ENERGY  DENSITY.    CHA0373 

DVOL=(XP-XI)*(XPx*2+XPXXI+XI**2)/3.D0  CHA0374 

RO   (I)=RONORM*(AROI  -  AROIP)/DVOL  CHA0375 

TENA(I)=RUNORM*(AROUI-AROUIP)/DVOL  CHA0376 

E    (I)  =  ENORM  *(AEI   -   AEIPVDVOL  CHA0377 

GO  TO  106  CHA0378 

105  CONTINUE  CHA0379 
C   UNIFORM  FLOW  REGION  CHA0380 

RO   (I)=RONORM*ROLIM  CHA0381 

TENA(I)=0.  CHA0382 

E    (I)=ENORM  x  ELIM  CHA0383 

106  CONTINUE  CHA0384 
U(I)=TENA(I)/RO(I)  CHA0385 
P(I)=(GAMA-1.DO)X(E(I)-0.5DOXRO(I)XU(I)XX2)  CHA0386 
PRINT  lll,I,CSOF,U(I),P(I),RO(I),E(I)  CHA0387 

111   FORMATC/1X, 'I,CS0F,U,P,R0,E=M4,L3,4D14.4)  CHA0388 

IF(CSOF)  GO  TO  109  CHA0389 

100  CONTINUE  CHA0390 

109  CONTINUE  CHA0391 
DO  4  1=2, LL  CHA0392 
DXSI(I)=(X(I+l)-X(I))XRO(I)  CHA0393 

<\           CONTINUE  CHA0394 

RETURN  CHA0395 

END CHAQ396 

SUBROUTINE  TAYLOR(GAMA)                          TA^  I    Oft  CHA0397 

IMPLICIT  REALX8(A-H,0-Z,$)                          J^^^  CHA0398 

C  CHA0399 

C   TAYLOR  --  SELF  SIMILAR  SPHERICAL  DETONATION  (CJ)  FLOW  FIELD  CHA0400 

C  CHA0401 

COMMON  /GGGG/G,G1,G2,G3,G<4,G5,G6,G7,G8,G9,G10  CHA0402 

COMMON  /PAR/RHO0,Q0,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,C0,U0  CHA0403 

COMMON  /GITN/NPO  CHA0404 

COMMON/GIT/ROLIM,ELIM,XGIT(200),ROGIT(200),ROUGIT(200),EGIT(200)   CHA0405 

CXXXXX5(«^5<5«^5(XXX^^XX^X^^^^^5(XXXXXXXXXXXXXXXXXX3€XXXXXXXXXX5(5(5(XXXK^XXXXXXXCHA0^06 

G=GAMA  CHA0407 

PRINT  101  CHA0408 

101  FORMATC'l1)  CHA0409 
PRINT  110  CHA0410 

110  FORMATUX, 'G.  I.  TAYLOR  SOLUTION.  N, PSI , U, C, X/AM, AT, AE= '// )  CHA0411 
CALL  INIDAT  CHA0412 
X  =  1.D0  CHA0413 
Y=0.  CHA0414 
U=U0  CHA0415 
C=C0  CHA0416 
AM=0.  CHA0417 
AT  =  0.  CHA0418 
AE=0.  CHA0419 
PSI=-DLOG(U)  CHA0420 
DO  1  N=1,NP0  CHA0421 
XGIT  (N)=X  CHA0422 
ROGIT  (N)=AM  CHA0423 
ROUGIT(N)=AT  CHA0424 
EGIT  (N)=AE  CHA0425 
PRINT  11,  N,PSI,U,C,X,AM,AT,AE  CHA0426 

11    FORMATUX, 14, 4D14.5/5X,3D14. 5)  CHA0427 

CALL  RUNGE(N,PSI,X,C,AM,AT,AE,PSIN,XN,CN,AMN,ATN,AEN)  CHA0428 

PSI=PSIN  CHA0429 

U=DEXP(-PSI)  CHA0430 

X=XN  CHA0431 

C=CN  CHA0432 

AM=AMN  CHA0433 

AT=ATN  CHA0434 

AE=AEN  CHA0435 

1     CONTINUE  CHA0436 

ROLIM=(C/C0)x*G3  CHA0437 

ELIM=G5x(C/C0)x*G4  CHA0438 

AM0=AM+(C/C0)XXG3XXXX3/3.D0  CHA0439 
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AM0=AM0*3.D0*(G+1.D0)/G  CHA0440 

AE0=AE+(G5X(C/C0)**G4+0.5D0*(C/C0)**G3*U*X2)XXXX3/3.D0  CHA0441 

AE0=AE0X6 .D0X(G+1.D0)X(GXX2-1.D0)/G  CHA0442 

PRINT  22,AMO,AEO  CHA0443 

!2    FORMATC///1X, 'MASS  AND  ENERGY  CHECK  (SHOULD  BE  1 . ) •//  CHA0444 

1  1X,,M0=',D17.8,5X,'E0=,,D17.8//)  CHA0445 

RETURN  CHA0446 

END  CHA0447 

SUBROUTINE   INlDAT  IHlDW  CHA0448- 

IMPLICIT  REALX8(A-H,0-Z,$)  CHA0449 

COMMON  /GGGG/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10  CHA0450 

COMMON  /PAR/RHO0,Q0,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,C0,U0  CHA0451 

COMMON  /GITN/NPO  CHA0452 

[xxxxxxkxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxcha  0^+53 

NP0  =  200  CHA0454 

PSIMAX=10.D0  CHA0455 

U0=1.D0/(G+1  DO)  CHA0456 

C0=1.D0-U0  CHA0457 

DPSI  =  PSIMAX/DFLOAT(NPO)  CHA0458 

G1=G-1.D0  CHA0459 

G2=G1/2.D0  CHA0460 

G3  =  2.D0/(G-1.D0)  CHA0461 

G4=2.D0XG/(G-1.D0)  CHA0462 

G5=G/((G+1.D0)XX2X(G-1.D0))  CHA0463 

RETURN  -iltr  c  CHA0464 

END KVri<*t.  rHAna^R 

SUBROUTINE  RUNGECN, PSI , X, C, AM, AT, AE, PSIN, XN, CN, AMN, ATN, AEN)  CHA0466 

IMPLICIT  REALX8(A-H,0-Z,$)  CHA0467 

COMMON  /GGGG/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10  CHA0468 

COMMON  /PAR/RHO0,Q0,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,C0,U0  CHA0469 

COMMON  /GITN/NPO  CHA0470 
;xxxxxxxxxxxxxxxxxxxxxkxkxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxchA0'+71 

H  =  DPSI  CHA0472 

H2=H/2.DO  CHA0473 

H6  =  H/6.D0  CHA0474 

CALL  DERIV(PSI,X,C,AM,AT,AE,  CHA0475 

1  DXDP1,DCDP1,DMDP1,DTDP1,DEDP1)  CHA0476 

CALL  DERIV(PSI+H2,X+H2*DXDP1,C+H2*DCDP1,AM,AT,AE,  CHA0477 

1  DXDP2,DCDP2,DMDP2,DTDP2,DEDP2)  CHA047S 

CALL  DERIV(PSI+H2,X+H2*DXDP2,C+H2*DCDP2,AM,AT,AE,  CHA0479 

1  DXDP3,DCDP3,DMDP3,DTDP3,DEDP3)  CHA0430 

CALL  DERIV(PSI+H,X+HXDXDP3,C+H*DCDP3,AM,AT,AE,  CHA04S1 

1  DXDP4,DCDP4,DMDP4,DTDP4,DEDP4)  CHA0432 

PSIN=PSI+H  CHA04S3 

XN=X+H6X(DXDP1+2.D0X(DXDP2+DXDP3)+DXDP4)  CHA0  484 

CN=C+H6*(DCDP1+2.D0*(DCDP2+DCDP3)+DCDP4)  CHA048  5 

AMN=AM+H6*(DMDP1+2.D0*(DMDP2+DMDP3)+DMDP4)  CHA0  486 

ATN=AT+H6*(DTDP1+2.D0*(DTDP2+DTDP3)+DTDP4)  CHA04S7 

AEN  =  AE+H6X(DEDP1  +  2.D0X(DEDP2  +  DEDP3)  +  DEDP<4)  CHA0483 
RETURN                                                 Tscnw,        CHA04S9 

END J/gKiy    CHA049Q 

SUBROUTINE  DERI V(  PSI  ,X,C,  AM,  AT,  AE,  DXDP,  DCDP,  DMDP,  DTDP,  DEDP)  CHA0491 

IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0492 

COMMON  /GGGG/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10  CHA0493 

COMMON  /PAR/RHO0,Q0,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,C0,U0  CHA0494 

COMMON  /GITN/NPO  CHA0495 
;xxxxxxxxxxxx^^¥xxxxxxxxxa«xxxxxxxxxxxxxxxxxxx^xxxxxxxxxxxxxxxxxxxxxxxxxcHA0'496 

U  =  DEXP(-PSI)  CHA0497 

DXDP=0.5D0*Xx(C-U+X)x(C+U-X)/Cx*2  CHA0  498 

DCDP=-G2XUX(X-U)/C  CHA0499 

DMDP=-(C/C0)**G3xX*x2*DXDP  CHA050  0 

DTDP=DMDP*U  CHA0501 

DEDP  =  -(G5*(C/C0)**G4+0.5D0*(C/C0)*XG3XUXX2;)*XXX2*DXDP  CHA0  50  2 

RETURN  CHA0503 

END . . CHA0504 

DOUBLE  PRECISION  FUNCTION  RATIO(X)  £/VTlO  CHA0505 

IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0506 

;xxxxxxxxxxxxxxxxxxxxx^x^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx5(xCHA0507 

RATIO=0.  CHA0503 

IF(X.LE.1.D-3)RETURN  CHA0509 

RATIO=2.D0/X  CHA0510 

RETURN  CHA0511 
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EN_D - CHAQ512 

DOUBLE  PRECISION  FUNCTION  CROSS(X)               CROSS  CHA0513 

IMPLICIT  REAL*8CA-H,0-Z,$)  CHA0514 

CXXXXKXXXXXXKXX5CXXXXXXXXXXXXXX5CXXXXXXXXXXXXXKXKXXXXXXXXXXXXXXXXXXXXXXXXXCHA0515 

C     CROSS=l.DO  CHA0516 

CR0SS=X*X2  CHA0517 

RETURN  CHA0518 

END . ("HAnsiQ 

SUBROUTINE  CYCEUL  CMC£-^L_     CHA0520 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0521 

2  US,PS,UIDOT,PIDOT,  CHA0522 
X             FIMZ,ZMDOT,  CHA0523 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0524 
IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0525 
DIMENSION  X(L),UCL),PCL),ROCL),GCL),ECL),DUCL),DPCL),DROCL),  CHA0526 

1  DG(L),DXSI(L),MIN(L),  CHA0527 

2  USCL),PSCL),UIDOTCL),PIDOTCL)  CHA0528 

3  ,TENACL),FIROCL),FIMCL),FIECL)  CHA0529 

4  ,GIPCL),VOLCL),ZCL),DZCL)  CHA0530 

5  ,FIMZCL),ZMDOTCL)  CHA0531 
COMMON  /AB/AC50)  CHA0532 
EQUIVALENCE  ( LL, A(2) ) , CT, AC3) ) , CDT, AC4) ) , (COLELA, A(13) )  CHA0533 
EQUIVALENCE  (KEYEK, A( 16 ) )  CHA0534 
EQUIVALENCE  (STAB, AC  18) ) , ( DTBA, AC  19) ) , ( DTKOD, A( 20) ) , CKDT, AC 21) )  CHA0535 
COMMON  /GAM/GAMA , NG, MU2 , Gl , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , Gl 0 , Gl 1  CHA0536 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA0537 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0538 
REAL*8  NG,MU2  CHA0539 
COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  CHA0540 
COMMON  /AZOV/ISAFA,NORIMN,USAF,PSAF,ROSAF,GSAF,ESAF,DPSAF  CHA0541 

1            ,DXSIL,DXSIR  CHA0542 

LOGICAL  NORIMN  CHA0543 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  CHA0544 

1  RSTARL,RSTARR,GSTARL,GSTARR,  CHA0545 

2  CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UWC6)  CHA0546 

3  ,LAMDAL,LAMDAR,RATEL, RATER, TEMPL , TEMPR, TEMPSL , TEMPSR  CHA0547 

4  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  CHA0543 
REALX8  LAMDAL,LAMDAR  CHA0549 
LOGICAL  HELEML,HELEMR  CHA0550 
COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR  CHA0551 

2  ,ASTARL,ASTARR,LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR  CHA0552 

3  ,RAT,SH  CHA0553 

4  ,BETACL,BETACR,DSDASL,DSDASR,DZDASL,DZDASR  CHA0554 
REALX8  LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR  CHA0555 
COMMON  /GRADS/DUDXIL,DPDXIL,DGDXIL,DRDXIL,DZDXIL,DSDXIL,  CHA0556 

1              DUDXIR,DPDXIR,DGDXIR,DRDXIR,DZDXIR,DSDXIR  CHA0557 

COMMON  /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN  CHA0558 

1  ,GIH  CHA0559 

2  ,FIH4,ZMD0TL,ZMD0TR  CHA0560 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET,  CHA0561 

1  RATE,TEMPC  CHA0562 
COMMON/PULS/PRESSC10),PULSE1C10),PULSE2C10),PULSE3C10),PULSE4C10)  CHA056  3 

DATA  ERRP/O.DO/  CHA0564 

DATA  NERRP/O/  CHA0565 

C     DATA  K0TZ/7777777777B/  CHA0566 
Cxxxxxxx*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCHA0  567 

DT2=DT/2.D0  CHA0568 

UXN=0.  CHA0573 

PXN=0.  CHA0574 

R0XN=0.  CHA0575 

ZXN=0.  CHA0576 

DO  1  1=2, L  CHA0577 

IM=I-1  CHA0578 

UXNM=UXN  CHA0579 

PXNM=PXN  CHA0580 

ROXNM=ROXN  CHA0581 

ZXNM=ZXN  CHA0582 

UL=UCIM)+0.5D0*DU(IM)  CHA0583 

PL=P(IM)+0.5D0*DP(IM)  CHA0584 

ROL=RO(IM)+0.5D0*DRO(IM)  CHA0585 

GL=DSQRTCGAMA*PLxROL)  CHA0586 

CL=GL/ROL  CHA0587 
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ZL=ZCIM)+0.5D0*DZ(IM) 

IF(ZL.GT.l.DO)  ZL=l.DO 

IFCZL.LT.O.   )  ZL=0. 

SL=PL/(G1XR0LXXGAMA) 

TEMPL=PL/ROL 

UR=U(I)-0.5D0XDUCI) 

PR=P(I)-0.5DOXDP(I) 

ROR=ROCI)-0.5D0XDRO(I) 

GR=DSQRT(GAMAXpRXROR) 

CR=GR/ROR 

ZR=ZCI)-0.5D0XDZ(I) 

IFCZR.GT.l.DO)  ZR=1.D0 

IFCZR.LT.O.   )  ZR=0. 

SR=PR/(G1XR0RXXGAMA) 

TEMPR=PR/ROR 

CALL  RIEMANCL,I,MIN) 

DUDXIL=DU(IM)/DXSI(IM) 

DPDXIL=DP(IM)/DXSI(IM) 

DRDXIL=DRO(IM)/DXSI(IM) 

DGDXIL=0.5D0XGLX(DPDXIL/PL+DRDXIL/ROL) 

DZDXIL=DZ(IM)/DXSI(IM) 

DSDXIL=SLX(DPDXIL/PL-GAMAXDRDXIL/ROL) 

DUDXIR=DU(I)/DXSI(I) 

DPDXIR=DP(I)/DXSI(I) 

DRDXIR=DRO(I)/DXSI(I) 

DGDXIR=0.5D0*GR*(DPDXIR/PR+DRDXIR/R0R) 
DZDXIR=DZ(I)/DXSI(I) 

DSDXIR=SRX(DPDXIR/PR-GAMAXDRDXIR/ROR) 
SH  =  CROSS(X(D) 
RAT  =  RATIO(X(D) 

CALL  MAGA(L,I,MIN) 


US(I)=USTAR 

PS(I)=PSTAR 

UIDOT(I)=DUIDT 

PIDOT(I)=DPIDT 

CALL  FLUXE(L,I,MIN) 

FIRO(I)=FIHl 
FIM  (I)=FIH2 
FIE  CI)=FIH3 
FIMZCI)=FIH4 
GIP(I)=GIH 
DU(IM)=UXN-UXNM 
DP(IM)=PXN-PXNM 
DRO(IM)=ROXN-ROXNM 
DZ(IM)=ZXN-ZXNM 
STATIONS  OUTPUT 

IF((I-42)X(I-6  2)X(I-82)X(I-102) 
NPU  =  0 

42)  NPU  =  1 

62)  NPU=2 

82)  NPU=3 
.102)NPU=4 
EQ.O)  CALL  SOF( 'FLUXE  90 


NE.O)  GO  TO  1 


IFCI.EQ 
IFCI.EQ 
IFCI.EQ 
IFCI.EQ, 

IF(NPU.EQ.O)  CALL  SOFC'FLUXE  90.   NPU.EQ.O') 

PRESS(NPU)=GIH+FIH2 

PULSE1CNPU)=PULSE1CNPU)+DT*GIH 

PULSE2CNPU)=PULSE2CNPU)+DT*CGIH+FIH2) 

PULSE3(NPU)=PULSE3(NPU)  +  DTXFIH1XCR0S3(X(D) 

PULSE4(NPU)=PULSE4(NPU)+DTXFIH2XCR05S(X(I)) 

CONTINUE 

AMTOT=0. 

ETOT=0. 

EKTOT=0. 

EPTOT=0. 

TENTOT=0. 

FI1=FIR0C2) 


CHA0588 
CHA0589 
CHA0590 
CHA0591 
CHA0592 
CHA0593 
CHA0594 
CHA0595 
CHA0596 
CHA0597 
CHA0598 
CHA0599 
CHA0600 
CHA0601 
CHA0602 
CHA0603 
CHA0604 
CHA0605 
CHA0606 
CHA0607 
CHA0608 
CHA0609 
CHA0610 
CHA0611 
CHA0612 
CHA0613 
CHA0614 
CHA0615 
CHA0616 
CHA0617 
CHA0618 
CHA0619 
CHA0620 
CHA0621 
CHA0622 
CHA0623 
CHA0624 
CHA0625 
CHA0626 
CHA0627 
CHA0623 
CHA0629 
CHA0630 
CHA0631 
CHA0632 
CHA0633 
CHA063<4 
CHA0635 
CHA0636 
CHA0637 
CHA0633 
CHA0639 
CHA0640 
CHA0641 
CHA0642 
CHA0643 
CHA0644 
CHA0645 
CHA0646 
CHA0647 
CHA0643 
CHA0649 
CHA0650 
CHA0651 
CHA0652 
CHA0653 
CHA0654 
CHA0655 
CHA0656 
CHA0657 
CHA0658 
CHA0659 
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FI2=FIM  (2)  CHA0660 

FI3=FIE  (2)  CHA0661 

FI4=FIMZ(2)  CHA0662 

GI2=GIP(2)  CHA0663 

SH=CR0SS(X(2))  CHA0664 

DO  2  1=2, LL  CHA0665 

IP=I+1  CHA0666 

FIM1=FI1  CHA0667 

FIM2=FI2  CHA0668 

FIM3=FI3  CHA0669 

FIM4=FI4  CHA0670 

GIM2=GI2  CHA0671 

SHM=SH  CHA0672 

FIl=FIRO(IP)  CHA0673 

FI2=FIM  (IP)  CHA0674 

FI3=FIE  (IP)  CHA0675 

FI4=FIMZ(IP)  CHA0676 

GI2=GIP  (IP)  CHA0677 

SH=CROSS(X(IP))  CHA0678 

DVOL=VOL(I)  CHA0679 

ROOLD=RO(I)  CHA0680 

POLD=P(I)  CHA0681 

EOLD=E(I)  CHA0682 

UOLD=U(I)  CHA0683 

ZOLD=Z(I)  CHA0684 

ZKODM=ZOLD*ROOLD  CHA0685 

TOLD=POLD/ROOLD  CHA0686 

DX=X(IP)-X(I)  CHA0687 

DTVOL=DT/DVOL  CHA0688 

C  CHA0689 

RO(I)=RO(I)-DTVOL*(SH*FIl-SHM*FIMl)  CHA06  90 

TENA(I)=TENA(I)-DTV0L*(SH*FI2-SHM*FIM2)-(DT/DX)X(GI2-GIM2)  CHA06  91 

E(I)=E(I)-DTV0L*(SH*FI3-SHM*FIM3)  CHA06  92 

U(I)=TENA(I)/RO(I)  CHA0693 

Z(I)=(ZK0DM-DTV0L*(SH*FI4-SHMXFIM4))/R0(I)  CHA06  94 

IF(Z(I).GT.1.D0)  Z(I)=1.D0  CHA0695 

IF(Z(I).LT.O.)  Z(I)=0.  CHA0696 

C  CHA0697 

UAV=U(I)  CHA0698 

ROAV=RO(I)  CHA0699 

EP=E(I)-0.5D0*ROAV*UAVX*2  CHA07  00 

IFCEP.GT.O.)  GO  TO  291  CHA0701 

NERRP=NERRP+1  CHA0702 

ERRP=ERRP+(l.D-3-EP)*DVOL  CHA07  03 

IF(ERRP.GT.0.24DO)  GO  TO  291  CHA0704 

EP=l.D-8  CHA0705 

291   CONTINUE  CHA0706 

IF(EP.LE.O.)  GO  TO  7001  CHA0707 

P(I)=G1*EP  CHA0708 

G(I)=DSQRT(GAMA*P(I)*RO(I))  CHA07  09 

C  CHA0710 

UPC=DABS(U(I))+G(I)/RO(I)  CHA0711 

DTI=STAB*DX/UPC  CHA0712 

IF(DTI.GT.DTBA)  GO  TO  29  CHA0713 

DTBA=DTI  CHA0714 

KDT=I  CHA0715 

29    CONTINUE  CHA0716 

DXSI(I)=RO(I)*DX  CHA0717 

ETOT=ETOT+E(I)*DVOL  CHA0718 

EPTOT=EPTOT+EP*DVOL  CHA0719 

AMTOT=AMTOT+RO(I)*DVOL  CHA0720 

TENTOT=TENTOT+TENA(I)*DVOL  CHA0721 

2     CONTINUE  CHA0722 

EKTOT=ETOT-EPTOT  CHA0723 

C  CHA0724 

IF(COLELA.EQ.0.)  GO  TO  200  CHA0725 

CALL  DCOLE(L,X,U  , DU  ,MIN,1)  CHA0726 

CALL  DC0LE(L,X,P,DP,MIN,2)  CHA0727 

CALL  DC0LE(L,X,R0,DR0,MIN,3)  CHA0728 

CALL  DC0LE(L,X,Z,DZ,MIN,4)  CHA0729 

200   CONTINUE  CHA0730 

CALL  BD0K1(L,X,U  ,  DU  ,MIN,1)  CHA0731 
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CALL  BD0K1(L,X,P,DP,MIN,2)  CHA0732 

CALL  BD0K1(L,X,R0,DR0,MIN,3)  CHA0733 

CALL  BD0K1CL,X,Z,DZ,MIN,4)  CHA0734 

PRINT  901,(NN,PRESS(NN),PULSE1(NN),PULSE2(NN),  CHA0735 

1  PULSE3(NN)/AMT0T,PULSE4(NN)/TENT0T,NN=1,4)  CHA0736 

01   F0RMAT(1X,2(»/M3,5D11.3,  »/')/)  CHA0737 

IF(DABS(T-A(5)).LT.l.D-6)  PRINT  911 , NERRP, ERRP  CHA0738 

11   F0RMATC//1X,  'NERRP,ERRP=M5,D15.5/)  CHA0739 

RETURN  CHA0740 

001  CONTINUE  CHA0741 
PRINT  7101,  I,ROAV,UAV,DROCI),DU(I),E(I),EP,ZNEW,ZNEW-1.DO,EPI     CHA0742 

101  F0RMATC//1X, "FROM  CYCEUL .   NEGATIVE  EP .   IN  CELL   I=f,I6//  CHA0743 

1  IX,  ,R0AV,UAV,DR0(I),DUCI)  =  ,,<tD18.8//  CHA0744 

2  1X,,E(I),EP,ZNEW,ZNEW-1,EPI=,,5D14.6//)  CHA0745 
CALL  PRINT  CHA0746 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0747 

2  US,PS,UIDOT,PIDOT,  CHA0748 
X             FIMZ,ZMDOT,  CHA0749 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0750 
CALL  SOFC'CYCEUL  7001,  NEGATIVE  EP  •  )  CHA0751 
RETURN  CHA0752 
END CHA0753 

SUBROUTINE  SAFAE SAfAF CHA0754 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,      n        ~  CHA0755 

2  US,PS,UIDOT,PIDOT,  CHA0756 
X             FIMZ,ZMDOT,  CHA0757 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0758 
IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0759 
DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  CHA076  0 

1  DG(L),DXSI(L),MIN(L),  CHA0761 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  CHA0762 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  CHA0763 

4  ,GIP(L),VOL(L),Z(L),DZ(L)  CHA0764 

5  ,FIMZ(L),ZMDOT(L)  CHA0765 
COMMON  /AB/AC50)  CHA0766 
EQUIVALENCE  ( LL , A( 2) ) , (T, A( 3) ) , ( DT, A(4 ) ) , (NCYC, AC  12) )  CHA0767 
EQUIVALENCE  (UGAL,A(15))  CHA0768 
COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  CHA076  9 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA0770 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0771 
REAL*8  NG,MU2  CHA0772 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET,  CHA077  3 

1  RATE,TEMPC  CHA0774 

C0MM0N/DIFFUS/U2,P2,R02,ARW  CHA077  5 

;XXX^^^XX^^^^X^X^^^^^^^X^X^5(^^XXX5(X5(X^X)(XX3(KXX5(3(X^X5«^3tXXX5<^X5(5«5(^X5(^^3«5(3<CHA0776 

CHA0777 

RIGID  B.C.  AT  1=2  CHA0778 

CHA0779 

U(l)=-U(2)  CHA0780 

P(1)=P(2)  CHA0731 

G(1)=G(2)  CHA0782 

R0(1)=R0(2)  CHA0733 

Z(1)=Z(2)  CHA0784 

DU(1)=DU(2)  CHA0785 

DP(1)=-DP(2)  CHA0786 

DG(1)=-DG(2)  CHA0787 

DR0(1)=-DR0(2)  CHA0738 

DXSI(1)=DXSI(2)  CHA0789 

CHA0790 

OUTFLOW  B.C.  AT  I=L  CHA0791 

CHA0792 

U(L)=U(LL)+DU(LL)/2.D0  CHA0793 

P(L)=P(LL)+DP(LL)/2.D0  CHA0794 

RO(L)=RO(LL)+DRO(LL)/2.D0  CHA07  95 

G(L)=G(LL)+DG(LL)/2.D0  CHA0796 

Z(L)=Z(LL)+DZ(LL)/2.D0  CHA0797 

DU(L)=0.  CHA0798 

DP(L)=0.  CHA0799 

DG(L)=0.  CHA0800 

DRO(L)=0.  CHA0801 

DZ(L)=0.  CHA0802 

DXSI(L)=DXSI(LL)  CHA0803 
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C  CHA0804 

RETURN  CHA0805 

END CHA0806 

SUBROUTINE  BDOK1  (  L  ,  X,  V,  DV,MIN,  NV)                RTjnkM  CHA0807 

IMPLICIT  REAL*8(A-H,0-Z,$)                      DJ^rx  CHA0808 

DIMENSION  X(L),V(L),DVCL),MIN(L)  CHA0809 

COMMON  /AB/AC50)  CHA0810 

EQUIVALENCE  ( LL , A(2) ) , ( KEYMON, A( 11 ) )  CHA0811 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX  CHA0812 

1             ,XMIN,XMAX,SMIN,SMAX,IVERSA  CHA0813 

COMMON  /M0NIT/CASEAV(4),NC14C4),NF16(6),  CHA0814 

1                NM0NU(4),NM0NP(4),NM0NR0(4),NM0NZC4)  CHA0815 

DIMENSION  NM0NV(4,4)  CHA0816 

EQUIVALENCE  ( NMONVC 1 , 1) , NMONUC 1) )  CHA0817 

DIMENSION  NAMEVC4)  CHA0818 

DATA  NAMEV/'U',  «P»,  'RO!,  'ZV  CHA0819 

DATA  EPS/1. D-9/  CHA0820 

CXXXXXXXXXXXXXXXXXXX3«XXXXXXXXXXXXXXXXXXXXXXXXXXX5«XXXXXXXXXXXXXXX3«XKXXXXXCHA0821 

GO  TO  (1,2,3,4),  NV  CHA0822 

1  AMIDA=(UMAX-UMIN)**2  CHA0823 
GO  TO  9  CHA0824 

2  AMIDA=(PMAX-PMIN)**2  CHA0825 
GO  TO  9  CHA0826 

3  AMIDA=(R0MAX-R0MIN)**2  CHA0827 
GO  TO  9  CHA0828 

4  AMIDA=1.D0  CHA0829 
GO  TO  9  CHA0830 

9     CONTINUE  CHA0831 

AMIDA=AMIDA*EPS**2  CHA0832 

EPSA=DSQRT(AMIDA)  CHA0833 

DO  29  1=2, LL  CHA0834 

ICAT=0  CHA0835 

IF(DABS(DV(I)).LE.EPSA)  DV(I)=0.  CHA0836 

IF(DVCI).EQ.O.)  GO  TO  29  CHA0837 

VLEFT=V(I)-0.5D0XDV(I)  CHA0838 

VRIGHT=V(I)+0.5D0*DV(I)  CHA0839 

VM=V(I-1)  CHA0840 

VP=V(I+1)  CHA0841 

SIGN=(VP-V(I))*(V(I)-VM)  CHA0842 

IF(SIGN.GT.-AMIDA)  GO  TO  22  CHA0843 

21  DV(I)=0.  CHA0844 
ICAT=1  CHA0845 
GO  TO  20  CHA0846 

22  CONTINUE  CHA0847 
SIGN=(VP-VM)*DV(I)  CHA0848 
IF(SIGN.GT.-AMIDA)  GO  TO  24  CHA0349 

23  DV(I)=0.5D0*(VP-VM)  CHA0850 
VLEFT=V(I)-0.5D0*DV(I)  CHA0851 
VRIGHT=V(I)+0.5D0*DV(I)  CHA0852 
ICAT=2  CHA0853 

24  SIGN=(VLEFT-VM)*DV(I)  CHA0854 
IF(SIGN.GT.-AMIDA)  GO  TO  26  CHA0855 

25  VLEFT=VM  CHA0856 
VRIGHT=2.D0*V(I)-VLEFT  CHA0857 
DV(I)=VRIGHT-VLEFT  CHA0858 
ICAT=3  CHA0859 

26  SIGN=(VP-VRIGHT)*DV(I)  CHA0860 
IF  (SIGN.GT.-AMIDA)  GO  TO  28  CHA0861 

27  VRIGHT=VP  CHA0862 
VLEFT=2.D0*V(I)-VRIGHT  CHA0863 
DV(I)=VRIGHT-VLEFT  CHA0864 
ICAT=3  CHA0865 

28  IF(DABS(DV(I)) ,LE.0.5D0*DABS(VP-VM))  GO  TO  31  CHA0866 

30  DV(I)=0.5D0*(VP-VM)  CHA0867 
ICAT=4  CHA0868 

31  CONTINUE  CHA0369 
20    CONTINUE  CHA0870 

IF  (DABS(DV(D).GT.EPSA)  GO  TO  40  CHA0871 

DV(I)=0.  CHA0872 

40    CONTINUE  CHA0873 

IF  (ICAT.GT.O)  NMONV(ICAT,NV)=NMONV(ICAT,NV)+l  CHA0874 

29  CONTINUE  CHA0877 
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RETURN  CHA0878 

END CHAflS79 

SUBROUTINE  DCOLEC  L  ,  X,  V,  DV,MIN,  NV)              T><LOLE  CHA0880 

IMPLICIT  REAL*8(A-H,0-Z,$)                    ^  CHA0881 

DIMENSION  X(L),V(L),DV(L),MIN(L)  CHA0882 

COMMON  /AB/AC50)  CHA0883 

EQUIVALENCE  (LL,A(2))  CHA0884 

XXXXXXX^XXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX3(XXXXXX3<XXXXXXXKXXXXXXXXCHA0885 

DO  1  1=2, LL  CHA0886 

IM=I-1  CHA0887 

IP=I+1  CHA0888 

DV(I)=0.5D0X(V(IP)-V(IM))  CHA0889 

1     CONTINUE  CHA0890 

RETURN  CHA0891 

END CHA0892 

SUBROUTINE  PRINT PQ\*\~r CHA0893 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,     r's"N  '     CHA0894 

2  US,PS,UIDOT,PIDOT,  CHA0895 
X             FIMZ,ZMDOT,  CHA0896 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  CHA0897 
IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0898 
DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  CHA0899 

1  DG(L),DXSI(L),MIN(L),  CHA0900 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  CHA0901 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  CHA0902 

4  ,GIP(L),VOL(L),Z(L),DZ(L)  CHA0903 

5  ,FIMZ(L),ZMDOT(L)  CHA0904 
COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  CHA0905 
COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  CHA0906 

1  RSTARL,RSTARR,GSTARL,GSTARR,  CHA0907 

2  CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UW(6)  CHA0908 

3  ,LAMDAL,LAMDAR,RATEL, RATER, TEMPL , TEMPR, TEMPSL , TEMPSR  CHA0909 

4  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  CHA0910 
REALX8  LAMDAL,LAMDAR  CHA0911 
LOGICAL  HELEML,HELEMR  CHA0912 
COMMON  /AB/AC50)  CHA0913 
EQUIVALENCE  ( LL , A(2) ) , (T, A( 3) ) , (NCYC, A( 12) ) , ( DT, A( 4) )  CHA0914 
EQUIVALENCE  (UGAL,A(15))  CHA0915 
C0MM0N/DIFFUS/U2,P2,R02,ARW  CHA0916 
COMMON/ DETO/QDET,PCJDET,RC J DET, UCJDET, DCJDET, PODET, ROODET,  CHA0917 

1            RATE,TEMPC  CHA0918 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  CHA0919 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,GZ1,G22,G23  CHA0920 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0921 
REAL*8  NG,MU2  CHA0922 
COMMON  /M0NIT/CASEAV(4),NC14(<4),NF16(6),  CHA0923 

1                NM0NU(4),NM0NP(4),NM0NR0C4),NM0NZ(4)  CHA0924 

DIMENSION  CASAVH4)  CHA0925 

LOGICAL  FULLPR  CHA0926 

XXXX^^^XXX^^X^5(XXXX^XX5«XXXXXXX^X3(XX^XXXXXXX^XX^XX^XXXXXX^^X^XXX)(XXX^^5(XCHA0  9  27 

FULLPR=.TRUE.  CHA0928 

PRINT  1  CHA0929 

1  FORMAT(lHl)  CHA0930 
PRINT  2,  T,DT,NCYC  CHA0931 

2  FORMAT(1X,10X, 'RESULTS  AT   T= ' , Dll . 5, 5X, ' DT= ' , Dl 1 . 5, 5X, ' NCYC= ' ,    CHA0932 
1                15//)  CHA0933 

PRINT  3,  AMTOT,ETOT,EKTOT,EPTOT,TENTOT  CHA0934 

3  FORMATC1X, ' AMTOT= • , D20 . 14, 2X, ' ETOT, EKTOT, EPTOT= • , 3D22 . 14/  CHA0955 
1        IX, «TENT0T=',D21.14//)  CHA0936 

4  FORMATUX,'   I','      X       ','      U       ','      P  ',      CHA0937 

1  '     RO       ','      G       »,'      Z  ',      CHA0938 

2  '     DU      ','     DP      ','    DRO      ',  CHA0939 

3  ■     DG      ','     DZ')  CHA0940 
44    FORMATC1X,'    «,»             »,'     US       ',»     PS  ',      CHA0941 

1  '    ZMDOT     ',•     FIMZ     ','   AMDOT  ',      CHA0942 

2  '    AMDOTN   ','    TEMP     ',»   ENTALP    ',  CHA0943 

3  •    AMACH    ','   ENTRO     ')  CHA0944 

5  FORMAT(IX)  CHA0945 
IF  (UGAL.NE.O.)  PRINT  6,  UGAL  CHA0946 

6  FORMATC/11X, 'INITIAL  VELOCITY  CORRESPONDS  TO  UGAL = ' , D15 . 6/ )  CHA0947 
DO  10  1=1, L  CHA0948 
IF  (MOD(I,10).NE.l)  GO  TO  11  CHA0949 
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PRINT  5  CHA0950 

PRINT  4  CHA0951 

PRINT  44  CHA0952 

PRINT  5  CHA0953 

11  CONTINUE  CHA0954 
PRINT  12,I,XCI),U(I),P(I),R0(I),GCI),Z(I),DUCI),DP(I),DR0(I),  CHA0955 

1           DG(I),DZ(I)  CHA0956 

12  F0RMAT(1X,I3,6D12.5,5D11.4)  CHA0957 
ENTRO=P(I)/RO(I)**GAMA  CHA0958 
IF(  .NOT.FULLPR)  GO  TO  131  CHA0959 
IF(I.EQ.l)  GO  TO  131  CHA0960 
IM=I-1  CHA0961 
UL=U(IM)+0.5*DU(IM)  CHA0962 
PL=P(IM)+0.5*DP(IM)  CHA0963 
ROL=R0(IM)+0.5*DRO(IM)  CHA0964 
GL=G(IM)+0.5*DG(IM)  CHA0965 
CL=GL/ROL  CHA0966 
ZL=Z(IM)+0.5*DZ(IM)  CHA0967 
IFCZL.LT.O.)  ZL=0.  CHA0968 
UR=U(I)-0.5*DU(I)  CHA0969 
PR=P(I)-0.5*DP(I)  CHA0970 
GR=G(I)-0.5*DGCI)  CHA0971 
ROR=RO(I)-0.5*DRO(I)  CHA0972 
CR=GR/ROR  CHA0973 
ZR=Z(I)-0.5*DZ(I)  CHA0974 
IFCZR.LT.O.)  ZR=0.  CHA0975 
IFCPL.LE.O.)  PL=l.D-8  CHA0976 
IFCPR.LE.O.)  PR=l.D-8  CHA0977 
CALL  RIEMAN(L,I,MIN)  CHA0978 
XI=X(I)  CHA0987 
RSTAR=RSTARL  CHA0988 
IF(USTAR.LT.O.)  RSTAR=RSTARR  CHA0989 
ZSTAR=ZL  CHA0990 
IFCUSTAR.LT.O.)  ZSTAR=ZR  CHA0991 
AMACH=USTAR/DSQRT(GAMA*PSTAR/RSTAR)  CHA0992 
AMDOT=RSTAR*USTAR*CROSSCXI)  CHA0993 
IFCI.NE.2)  GO  TO  132  CHA0994 
AMDOT0=AMDOT  CHA0995 
IF(DABSCAMDOTO) .LT.l.D-12)  AMDOTO=1.DO  CHA0996 

132   CONTINUE  CHA0997 

AMDOTN=AMDOT/AMDOT0  CHA0998 
ENTALP=(GAMA/(GAMA-l.DO))*PSTAR/RSTAR+0.5DO*USTAR**2+QDET*ZSTAR    CHA0999 

ARW=1.D0  CHAIOOO 

TEMP=PSTAR/(RSTAR*ARW)  CHA1001 

PRINT  13,US(I),PS(I),  CHA1002 
1         ZMDOT(I),FIMZ(I),AMDOT,AMDOTN,TEMP,ENTALP,AMACH,ENTRO     CHA1003 

13  F0RMAT(4X,12X,5D12.5,6D11.4)  CHA1004 
131  CONTINUE  CHA1005 
10    CONTINUE  CHA1006 

C   JOB  STATISTICS  CHA1007 

DO  40  1=1,4  CHA1008 

CASAV1(I)=0.  CHA1009 

IF  (NC14(I) .NE.O)  CASAV1(I)=CASEAV(I)/DFL0AT(NC14(I))  CHA1010 

40    CONTINUE  CHA1011 

PRINT  30  CHA1012 

30  F0RMAT(///1X,10( '*'),3X, 'JOB  STATISTICS ', 3X, 10( ' *' )// )  CHA1013 
PRINT  31,(NC14(I),I=1,4)  CHA1014 

31  FORMATUX, 'NO.  OF  VARIOUS  CASES  IN  RIEMAN  SOLVER  NC14C NCASE) = « ,  CHA1015 
1             4110)  CHA1016 

PRINT  301,  (CASAV1(I),I=1,4)  CHA1017 

301   F0RMATC/1X, 'AVERAGE  NUMBER  OF  ITERATIONS  IN  RIEMAN  SOLVER',  CHA1018 

1        IX,*   CASAV1(NCASE)=',4(F6.2,4X))  CHA1019 

PRINT  32, (NF16(I), 1=1,6)  CHA1020 

32  F0RMATC/1X, 'NO.  OF  VARIOUS  FLUX  CASES  NF16 (NFLUX) = ' , 6110 )  CHA1021 
ICAT0=4  CHA1022 
PRINT  33,(NMONU(I),I=l,ICAT0),(NMONP(I),I=l,ICAT0),  CHA1023 

1          (NMONRO(I),I=1,ICATO),(NMONZ(I),I=1,ICATO)  CHA1024 

33  F0RMATC/1X, 'NO.  OF  MONOTONICITY  INTERVENTIONS  FOR  EACH  VAR . ' ,  CHA1025 
1  IX, 'IN  EACH  CATEGORY.'/  CHA1026 
1  1X,'NM0NU  (ICAT)=',4I10/  CHA1027 
1  1X,'NM0NP  (ICAT)=',4I10/  CHA1028 
1        IX, 'NMONRO(ICAT)=',4I10/  CHA1029 
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RETURN 
END 


1X,'NM0NZ  (ICAT)=',4I10/) 


CHA1030 
CHA1031 
CHA1032 


SUBROUTINE  SOF(ISTOP) 

IMPLICIT  REAL*8(A-H,0-Z,$) 

DIMENSION  ISTOP(l) 

PRINT  1,  ISTOP 

FORMAT(//1X,3H**X,3X,20A4,3X,3HXXX///) 

PRINT  1 

XX=-1.D0 

YY=DSQRT(XX) 

STOP 

END      


SOf 


CHA1034 
CHA1035 
CHA1036 
CHA1037 
CHA1038 
CHA1039 
CHA1040 
CHA1041 
CHA1042 
CHAlfl^ 


R I E  MAN 


SUBROUTINE  RIEMAN( L , I , MIN) 

IMPLICIT  REAL*8(A-H,0-Z,$) 

DIMENSION  MIN(L) 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1  RSTARL,RSTARR,GSTARL,GSTARR, 

2  CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UW(6) 

3  , LAMDAL,LAMDAR,RATEL, RATER, TEMPL , TEMPR,TEMPSL ,TEMPSR 

4  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
REAL*8  LAMDAL,LAMDAR 
LOGICAL  HELEML,HELEMR 
COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR 

2  , ASTARL , ASTARR, LAMDSL , LAMDSR, DSDAL , DSDAR, DZDAL, DZDAR 

3  ,RAT,SH 

4  , BETACL , BETACR, DSDASL , DSDASR, DZDASL , DZDASR 
REALX8  LAMDSL , LAMDSR, DSDAL , DSDAR, DZDAL , DZDAR 
COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2  , G24 , G25 , G26 , G27 , G28 , G29 , G30 , G31 , G32, G33 , G34 ,  G35 
REAL*8  NG,MU2 
COMMON  /AB/A(50) 
COMMON  /M0NIT/CASEAV(4),NC14(4),NF16(6), 

1  NM0NU(4),NM0NP(4),NM0NR0(4),NM0NZ(<+) 

(XXXXXXXXXXXXxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCHA1334 
DATA  NMAX/63/  CHA1335 

DATA  EPS/1. D-8/  CHA1336 

DATA  NTRY/0/  CHA1337 

(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX^^X5«XCHA1338 
UW(6)=1.D20  CHA1339 

WL=0.  CHA1340 

WR=0.  CHA1341 

ZETAL=PL**G8  CHA1342 

ZETAR=PR**G8  CHA1343 

CLG=CL/GAMA  CHA1344 

CRG=CR/GAMA  CHA1345 

ZSTARL=ZL  CHA1346 

ZSTARR=ZR  CHA1347 

IF  (ZETAL.LT.ZETAR)  GO  TO  102  CHA1343 

LEFT  PRESSURE  IS  HIGHER  CHA1549 


CHA1310 
CHA1311 
CHA1312 
CHA1313 
CHA1314 
CHA1315 
CHA1316 
CHA1317 
CHA1318 
CHA1319 
CHA1320 
CHA1321 
CHA1322 
CHA1323 
CHA1324 
CHA1325 
CHA1326 
CHA1327 
CHA1328 
CHA1329 
CHA1330 
CHA1331 
CHA1332 
CHA1333 


)1 


D0+G6*EVERR) 


IF 
IF 
IF 
GO 
RIGHT 
32 


CONTINUE 

EVERR=CPL-PR)/PR 

USR=UR+CRG*EVERR/DSQRT(1 

SRR=USR 

UEL=UL-G7XCLX(ZETAR-ZETAL)/ZETAL 

SLL=UEL 

NL=2 

NR  =  2 

IF  (USR.GE.UL)  NL=1 

(UEL.LE.UR)  NR=1 

(DABS(EVERR) .LT.EPS)  GO  TO  100 

(NL.EQ.2.AND.NR.EQ.1)  GO  TO  7001 

TO  100 

PRESSURE  IS  HIGHER 
CONTINUE 
EVERL=(PR-PL)/PL 
USL=UL-CLG*EVERL/DSQRT(1 
SLL=USL 
UER=UR+G7*CR*(ZETAL-ZETAR)/ZETAR 


D0+G6XEVERL) 


CHA1350 
CHA1351 
CHA1352 
CHA1353 
CHA1354 
CHA1355 
CHA1356 
CHA1357 
CHA1353 
CHA1359 
CHA1360 
CHA1361 
CHA1362 
CHA1363 
CHA1364 
CHA1365 
CHA1366 
CHA1367 
CHA1368 
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100 


SRR=UER 

NL  =  2 

NR  =  2 

IF  (UER.GE.UL)  NL=1 

IF  (USL.LE.UR)  NR=1 

IF  (DABS(EVERL).LT.EPS)  GO  TO  100 

IF  (NL.EQ.1.AND.NR.EQ.2)  GO  TO  7001 

GO  TO  100 

CONTINUE 

IF  (NL.EQ 

IF  (NL.EQ 


AND.NR.EQ.2)  NCASE=1 

AND.NR.EQ.2)  NCASE=2 
IF  (NL.EQ.2.AND.NR.EQ.1)  NCASE=3 
IF  (NL.EQ.l.AND.NR.EQ.l)  NCASE=4 

IF(DABS(PL-PR)+DABS(UL-UR) .LT.EPSX(PMAX-UMIN))  NCASE=4 
UMIDA=EPS*DMAX1(CL,CR) 
DUDZL=-G7*CL/ZETAL 
DUDZR=  G7XCR/ZETAR 

ZETA=(-(UR-UL)+ZETAR*DUDZR-ZETAL*DUDZL)/(DUDZR-DUDZL) 
IF  (ZETA.LE.O.)  GO  TO  7002 
N  =  0 

GO  TO  (1,2,3,4),  NCASE 
:   THE  CASE   ES 

1  ITYPE=NCASE 
HELEML=. FALSE. 
HELEMR=.TRUE. 

11    N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 

ZETAF=ZETA 

UEL=UL-G7*CL*(ZETAF-ZETAL)/ZETAL 

PPR=(ZETAF/ZETAR)**NG 

EVERR=PPR-1.D0 

SQRR=DSQRT(1.D0+G6*EVERR) 

USR=UR+CRG*EVERR/SQRR 

DU=UEL-USR 

IF  (DABS(DU).LE.UMIDA)  GO  TO  10 

DUDZR=NG*CRG*(PPR/ZETAF)*( 1 . D0+G9*EVERR)/SQRR**3 

ZETA=ZETAF+DU/(DUDZR-DUDZL) 

GO  TO  11 
10    CONTINUE 

USTAR=(UEL+USR)/2.D0 

IF(DABS(USTAR) . LT . EPS*UMAX)    USTAR=0. 

PSTAR=PPR*PR 

CSTARL=CL+(UL-USTAR)/G7 

RSTARL=GAMA*PSTAR/CSTARL**2 

GSTARL=CSTARL*RSTARL 
:   EQU.  NO.   69.01  OF  THE  BOOK  BY  COURANT-FRIEDRICHS . 

HWR=G11*(USTAR-UR)*R0R 

WR=WWR+DSQRT(GR**2+WHR**2) 

RSTARR=ROR*HR/(WR-ROR*(USTAR-UR)) 

GSTARR=DSQRT(GAMA*PSTAR*RSTARR) 

CSTARR=GSTARR/RSTARR 

WRE=WR/ROR+UR 

UH(1)=UL-CL 

UW(2)=USTAR-CSTARL 

UW(3)=USTAR 

UW(4)=WRE 

UH(5)=WRE 

GO  TO  5 
:   THE  CASE   SS 

2  ITYPE=NCASE 
HELEML=.TRUE. 
HELEMR=.TRUE. 

21    N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 

ZETAF=ZETA 

PF=ZETAF**NG 

PPL=PF/PL 

PPR=PF/PR 

EVERL=PPL-1 .DO 

EVERR=PPR-1.D0 

SQRL=DSQRT(1.D0+G6*EVERL) 

SQRR=DSQRT(1 . D0+G6XEVERR) 
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20 

D0+G9*EVERL)/SQRL**3 

D0+G9*EVERR)/SQRR**3 


USTAR=0 


USL=UL-CLG*EVERL/SQRL 
USR=UR+CRG*EVERR/SQRR 
DU=USL-USR 

IF  (DABS(DU).LE.UMIDA)  GO  TO 
DUDZL=-NG*CLG*(PPL/ZETAF)*(1 
DUDZR=  NG*CRG*(PPR/ZETAF)X(1 
ZETA=ZETAF+DU/(DUDZR-DUDZL) 
GO  TO  21 
I    CONTINUE 

USTAR=(USL+USR)/2.D0 
IF(DABS(USTAR) . LT . EPS*UMAX) 
PSTAR=(PPL*PL+PPR*PR)/2.D0 
WWR=G11*(USTAR-UR)*R0R 
WR=WWR+DSQRT(GR**2+WWR**2) 
WWL=-G11*(USTAR-UL)*R0L 
WL=WWL+DSQRT(GL**2+WWL**2) 
RSTARL  =  ROL*WL/(WL  +  ROL*(USTAR-UD) 

RSTARR=ROR*WR/(WR-ROR*(USTAR-UR)) 

GSTARL=DSQRTCGAMA*PSTAR*RSTARL) 

GSTARR=DSQRT(GAMA*PSTAR*RSTARR) 

CSTARL=GSTARL/RSTARL 

CSTARR=GSTARR/RSTARR 

WLE=-WL/ROL+UL 

WRE=WR/ROR+UR 

UW(1)=WLE 

UW(2)=WLE 

UW(3)=USTAR 

UW(4)=WRE 

UW(5)=WRE 

GO  TO  5 
THE  CASE   SE 

ITYPE=NCASE 

HELEML=.TRUE. 

HELEMR=. FALSE. 

N  =  N+1 

IF  (N.GT.NMAX)  GO  TO  7003 

ZETAF=ZETA 

UER=UR+G7*CR*(ZETAF-ZETAR)/ZETAR 

PPL=(ZETAF/ZETAL)**NG 

EVERL=PPL-1.D0 

SQRL=DSQRT(1 .  D0  +  G6*EVERL ) 

USL=UL-CLG*EVERL/SQRL 

DU=USL-UER 

IF  (DABS(DU) .LE.UMIDA)  GO  TO  30 

DUDZL=-NG*CIG*(PPL/ZETAF)X(1.D0+G9*EVERL)/SQRLX*3 

ZETA=ZETAF+DU/(DUDZR-DUDZL) 

GO  TO  31 
I    CONTINUE 

USTAR=(USL+UER)/2.D0 

IFCDABS(USTAR) . LT . EPS*UMAX)  USTAR=0. 

PSTAR=PPL*PL 

CSTARR=CR-(UR-USTAR)/G7 

RSTARR=GAMA*PSTAR/CSTARR**2 

GSTARR=CSTARR*RSTARR 

WWL=-G11*(USTAR-UL)*R0L 

WL=WHL+DSQRT(GL*X2+WHLXX2) 

WLE=-WL/ROL+UL 

RSTARL  =  ROL*WL/(WL+ROL*(USTAR-UD) 

GSTARL=DSQRT(GAMA*PSTAR*RSTARL) 

CSTARL=GSTARL/RSTARL 

UW(1)=WLE 

UW(2)=WLE 

UW(3)=USTAR 

UW(4)=USTAR+CSTARR 

UH(5)=UR+CR 

GO  TO  5 
THE  CASE   EE 

ITYPE=NCASE 

HELEML=. FALSE. 

HELEMR=. FALSE. 

PSTAR=ZETA*XNG 

USTAR=UL-G7XCL*(ZETA-ZETAL)/ZETAL 
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IF(DABS(USTAR).LT.EPS*UMAX)  USTAR=0.  CHA1513 

CSTARL=CL  +  (UL-USTAR)/G7  CHA1514 

CSTARR=CR-(UR-USTAR)/G7  CHA1515 

RSTARL=GAMA*PSTAR/CSTARL**2  CHA1516 

RSTARR=GAMA*PSTAR/CSTARR**2  CHA1517 

GSTARL=RSTARL*CSTARL  CHA1518 

GSTARR=RSTARR*CSTARR  CHA1519 

UW(1)=UL-CL  CHA1520 

UW(2)=USTAR-CSTARL  CHA1521 

UW(3)=USTAR  CHA1522 

UW(4)=USTAR+CSTARR  CHA1523 

UW(5)=UR+CR  CHA1524 

N=l  CHA1525 

GO  TO  5  CHA1526 

5  CONTINUE  CHA1527 
DO  6  K=l,6  CHA1528 
NFLUX=K  CHA1529 
IF  (UW(K).GE.O.)  GO  TO  61  CHA1530 

6  CONTINUE  CHA1531 
NFLUX=6  CHA1532 

61    CONTINUE  CHA1533 

NC14(NCASE)=NC14(NCASE)+1  CHA1537 

CASEAV(NCASE)=CASEAV(NCASE)+DFLOAT(N)  CHA1538 

NF16(NFLUX)=NF16(NFLUX)+1  CHA1539 

IF(NTRY.GE.2)G0  TO  666  CHA1540 

IFCI.NE.2.AND.I.NE.L)  GO  TO  666  CHA1541 
PRINT  667,I,NFLUX,NCASE,PL,UL,R0L,PR,UR,R0R,USTAR,PSTAR,RSTARL,    CHA1542 

1  RSTARR,(KK,UW(KK),KK=1,6)  CHA1543 

667    FORMATC/1X,  '  I ,  NFLUX,  NCASE=  ' ,  3I5/1X,  "PL,  UL  ,  ROL,  PR,  UR,  ROR=  f  ,6D12  .  4/  CHA1544 

1  IX, 'USTAR,PSTAR,RSTARL,RSTARR=',4D13.4/  CHA1545 

2  lX,,KK,UW(KIO  =  ',6(I4,2X,D13.4)/)  CHA1546 
NTRY=NTRY+1  CHA1547 

666   CONTINUE  CHA1548 

RETURN  CHA1549 

7001  CONTINUE  CHA1550 
PRINT  7101,  PL,UL,PR,UR,ZETAL,ZETAR,SLL,SRR,NL,NR,I  CHA1551 

7101  F0RMAT(//1X, 'FROM  RIEMAN.   AN  IMPOSSIBLE  CASE  OF  EXPANSION/SHOCK'  CHA1552 

1  //1X, 'PL,UL,PR,UR=',4D25.14//  CHA1553 

2  IX, 'ZETAL,ZETAR,SLL,SRR=',4D25.14//  CHA1554 

3  IX, 'NL,NR,I=',3I10//)  CHA1555 
CALL  SOFC7001')  CHA1556 

7002  CONTINUE  CHA1557 
PRINT  7102,  ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR,N,NCASE,I     CHA1558 


7102  FORMATC//1X, 

1  IX, 

2  IX, 

3  IX, 

4  IX, 

5  IX, 


FROM  RIEMAN.   NEGATIVE  PRESSURE  AT  THE  INTERSECTION ', CHA1559 

OF  L  AND  R  EXPANSION  BRANCHES'//  CHA1560 

IT  MEANS  THAT  A  CAVITATION  TENDS  TO  FORM.   THIS',  CHA1561 

POSSIBILITY  IS  EXCLUDED  IN  PRESENT  VERSION'//  CHA1562 

ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR=',9D10.3//   CHA156  3 

N,NCASE,I=* ,3110//)  CHA1564 

CALL  SOF( '7002')  CHA1565 

7003  CONTINUE  CHA1566 

PRINT  7103,  I,N,NCASE,DU,UMIDA,EPS,PL,UL,PR,UR,  CHA1567 

1            ZETA,ZETAF,ZETAL,ZETAR,DUDZL,DUDZR  CHA1568 

7103  FORMATC//1X, 'FROM  RIEMAN.   NUMBER  OF  ITERATIONS  EXCEEDED.1//  CHA1569 

1  IX, 'I,N,NCASE,DU,UMIDA,EPS=',3I6,3D18.6//  CHA1570 

2  IX, 'PL,UL,PR,UR,ZETA,ZETAF=«,6D18.10//  CHA1571 

3  IX, 'ZETAL,ZETAR,DUDZL,DUDZR=',4D18.10//)  CHA1572 
CALL  SOFC7003')  CHA1573 
RETURN  CHA1574 
END  CHA1575 

C$OPTIONS  LIST CHA1576 

SUBROUTINE  MAGACL, I, MIN)                          M/L/1A-  CHA1577 

IMPLICIT  REAL*8(A-H,0-Z,$)                         '  CHA1578 

DIMENSION  MIN(L)  CHA1579 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  CHA158  0 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA1531 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA1582 
REALX8  NG,MU2  CHA1583 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET,  CHA1584 

1            RATE,TEMPC  CHA1585 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  CHA1586 

1              RSTARL,RSTARR,GSTARL,GSTARR,  CHA1587 
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2  CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UW(6)  CHA1588 

3  ,LAMDAL,LAMDAR,RATEL, RATER, TEMPL , TEMPR, TEMPSL , TEMPSR  CHA1589 

4  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  CHA1590 
REALX8  LAMDAULAMDAR  CHA1591 
LOGICAL  HELEML,HELEMR  CHA1592 
COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR  CHA1593 

2  ,ASTARL,ASTARR,LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR  CHA1594 

3  ,RAT,SH  CHA1595 

4  ,BETACL,BETACR,DSDASL,DSDASR,DZDASL,DZDASR  CHA1596 
REAL*8  LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR  CHA1597 
COMMON  /GRADS/DUDXIL,DPDXIL,DGDXIL,DRDXIL,DZDXIL,DSDXIL,  CHA1598 

1              DUDXIR,DPDXIR,DGDXIR,DRDXIR,DZDXIR,DSDXIR  CHA1599 

COMMON  /AB/AC50)  CHA1600 

REALX8  LU,LP,LRO,LLAMDA  CHA1601 

DATA  EPS/1. D-6/  CHA1602 

XXXXXXX3«XXXXXXXXXX5(X5(XXXXX5<XKKXXXXXXXXX3(5(XXXX3(X3(XXXX5(XX5CX^)(XXX3(3(5(?(XXXXXCHA16  03 

WE  HERE  SOLVE  FOR  THE  TIME-DERIVATIVES  ALONG  THE  CONTACT  SURFACE,     CHA1604 

NAMELY   DUIDT,DPIDT.   FROM  THESE  WE  ALSO  OBTAIN  THE  OTHER  CHA1605 

TIME-DERIVATIVES  (SEE  COMMON  /STEP1/).  CHA1606 
WE  COMPUTE  THE  COEFFICIENTS  FOR  TWO  EQUATIONS  FOR  DUIDT,DPIDT.   THESECHA1607 

ARE         AAL*DUIDT+BBL*DPIDT=DDL  CHA1608 

AAR*DUIDT+BBR*DPIDT=DDR  CHA1609 

XKXXXXXXXXXXX3(X3(XXXXXX)(XXXXXXXXXXXXXXXXXXXXXXXXXX3<XXXXXXXXXXXXXXXX3O(XX5(CHA1610 

IF(SH.LE.EPS)RAT=0.  CHA1611 

CHA1612 

LEFT  SIDE  OF  CONTACT  CHA1613 

CHA1614 

IF  (.NOT.HELEML)  GO  TO  12  CHA1615 

11  CONTINUE  CHA1616 
LEFT  SHOCK  CHA1617 

DP=PSTAR-PL  CHA1618 

DU=USTAR-UL  CHA1619 

Z2=0.5D0/(PSTAR+MU2*PL)  CHA1620 

LU=DU*(0 .5DOXROL+MU2*Z2*GLX*2)-GL*X2/WL-WL  CHA1621 

LRO=-0.5D0*DP/ROL  CHA1622 

LP=-2.D0-MU2*Z2*DP  CHA1623 

AAL=2.D0-Z2*DP  CHA1624 

BBL=Z2*DU+HL/GSTARL**2+1 . DO/WL  CHA1625 

DDL=LUXDUDXIL+LRO*DRDXIL+LP*DPDXIL  CHA16  26 

DDL=DDL-WL*USTAR*RAT/RSTARL  CHA16  27 

1   +UL*RAT*(-GAMA*PL/WL+DU*(GAMA*PLXMU2*Z2+0.5D0))  CHA1628 

GO  TO  10  CHA1629 

12  CONTINUE  CHA1630 
LEFT  RAREFACTION  CHA1631 

A1=DUDXIL+DPDXIL/GL  CHA1632 

BETA=GSTARL/GL  CHA1633 

SQB=DSQRT(BETA)  CHA1634 

ASTARL=Al-(CL/(G15*SL))*DSDXIL*(BETAx*G5-l.D0)  CHA16  35 

AAL=1.D0  CHA1636 

BBL=1.D0/GSTARL  CHA1637 

DDL=-GSTARL*ASTARL/SQB  CHA1638 

DSDAL=DSDXIL  CHA1639 

DZDAL=DZDXIL  CHA1640 

DSDASL=DSDXIL*SQB  CHA1641 

DZDASL=DZDXIL*SQB  CHA1642 

GEOM  =  RAT*((GAMA-l.D0)*UL  +  2.D0XCD*  CHA16  43 

1  (BETA**G13-l.D0)/(ROL*(GAMA-3.D0))  CHA1644 

1   -4.D0*RAT*CL*(BETA**G14-1 . DO )/ ( ROL*( 3 . D0*GAMA-5 . DO ) )  CHA1645 

ASTARL=ASTARL-GEOM  CHA1646 

EVER1=  GSTARL*GEOM/SQB  CHA1647 

EVER2=-RAT*USTAR*CSTARL  CHA1643 

DDL=DDL+EVER1+EVER2  CHA1649 

GO  TO  10  CHA1650 

10    CONTINUE  CHA1651 

CHA1652 

RIGHT  SIDE  OF  CONTACT  CHA1653 

CHA1654 

IF  ( .NOT.HELEMR)  GO  TO  22  CHA1655 

21    CONTINUE  CHA1656 

RIGHT  SHOCK  CHA1657 

DP=PSTAR-PR  CHA1658 

DU=USTAR-UR  CHA1659 
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Z2=0.5D0/(PSTAR+MU2*PR) 

LU=DU*(0.5DO*ROR+MU2*Z2*GR**2)+GR*X2/WR+WR 
LR0=-0.5D0XDP/R0R 

LP=-2.D0-MU2XZ2XDP 

AAR=2.D0-Z2XDP 

BBR=Z2XDU-WR/GSTARRXX2-1 . DO/WR 

DDR=LUXDUDXIR+LROXDRDXIR+LPXDPDXIR 

DDR=DDR+WRXUSTARXRAT/RSTARR 

1    +URXRAT*(GAMAXPR/WR+DUX(GAMAXPRXMU2*Z2+0.5D0)) 
GO  TO  20 

22    CONTINUE 

C   RIGHT  RAREFACTION 

A1=DUDXIR-DPDXIR/GR 

BETA=GSTARR/GR 

SQB=DSQRTCBETA) 

ASTARR=A1+CCR/CG15XSR))XDSDXIRXCBETAXXG5-1.D0) 

AAR=1.D0 

BBR=-1.D0/GSTARR 

DDR=GSTARR*ASTARR/SQB 

DSDAR=DSDXIR 

DZDAR=DZDXIR 

DSDASR=DSDXIRXSQB 

DZDASR=DZDXIRXSQB 

GEOM=RAT*(-(GAMA-l.D0)*UR+2.D0*CR)X(BETAXXG13-l.D0) 

1  /CRORXCGAMA-3.D0)) 

2  -4 . D0XRAT*CRX(BETAXXG14-1 ' D0)/(RORX(3 . D0XGAMA-5 .  DO) ) 
ASTARR=ASTARR+GEOM 

EVER1=GSTARRXGE0M/SQB 
EVER2=RATXUSTARXCSTARR 
DDR=DDR+EVER1+EVER2 
GO  TO  20 
20    CONTINUE 

DET=AALXBBR-AARXBBL 

DUIDT=(DDLXBBR-DDRXBBL)/DET 

DPIDT=-(DDLXAAR-DDRXAAL)/DET 

DRIDTL=DPIDT/CSTARLXX2 

DRIDTR=DPIDT/CSTARRXX2 

RETURN 

END 


SUBROUTINE  FLUXEC L , I , MIN) 

IMPLICIT  REAL*8CA-H,0-Z,$) 

DIMENSION  MIN(L) 

COMMON  /AB/A(50) 

EQUIVALENCE  ( DT, AC  4) ) , ( NCYC, AC  12) ) 

COMMON  /GAM/GAMA , NG , MU2 , Gl , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , Gl 0 , Gl 1 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REALX8  NG,MU2 

COMMON  /GRADS/DUDXIL, DPDXIL , DGDXIL , DRDXIL , DZDXIL , DSDXIL, 
1  DUDXIR,DPDXIR,DGDXIR,DRDXIR,DZDXIR,DSDXIR 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1  RSTARL,RSTARR,GSTARL,GSTARR, 

2  CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UWC6) 

3  , LAMDAL,LAMDAR,RATEL, RATER, TEMPL , TEMPR, TEMPSL 

4  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
REALX8  LAMDAL,LAMDAR 

LOGICAL  HELEML,HELEMR 

COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR 

2  ,  ASTARL , ASTARR, LAMDSL , LAMDSR, DSDAL ,  DSDAR, DZDAL , DZDAR 

3  ,RAT,SH 

4  , BETACL , BETACR, DSDASL , DSDASR, DZDASL , DZDASR 
REAL 5(8  LAMDSL  ,  LAMDSR,  DSDAL  ,  DSDAR,  DZDAL  ,  DZDAR 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET, 

1  RATE,TEMPC 

COMMON  /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN 

1  ,GIH 

2  ,FIH4,ZMD0TL,ZMD0TR 
REALX8  LAMDAO 

CXXXXXXXXX*XX*XXXXXX**XXXXXXXXXXXX***X*X**XX*X*X**X*X*X***X*XX*X 

C   RO,U,P,Z  AND  THEIR  CXI,T)  DERIVATIVES  AT  EULERIAN  POINT  X=XCI 

CXXXX*X*XXXXXX*XX**X*X*XX**X*X*X*XX*XX*X*XXX*XXX****XXX*X*****X* 

DT2=DT/2.D0 
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CHA1665 
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CHA1673 
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CHA1675 
CHA1676 
CHA1677 
CHA1678 
CHA1679 
CHA1680 
CHA1681 
CHA1682 
CHA1683 
CHA1684 
CHA1685 
CHA1686 
CHA1687 
CHA1688 
CHA1689 
CHA1690 
CHA1691 
CHA1692 
CHA1693 
CHA1694 
CHA1695 
CHA1696 
CHA1697 
CHA1698 


FLUXt 


CHA1699 
CHA1700 
CHA1701 
CHA1702 
CHA1703 
CHA1704 
CHA1705 
CHA1706 
CHA1707 
CHA1708 
CHA1709 
CHA1710 
CHA1711 
CHA1712 
,TEMPSR  CHA1713 
CHA1714 
CHA1715 
CHA1716 
CHA1717 
CHA1718 
CHA1719 
CHA1720 
CHA1721 
CHA1722 
CHA1723 
CHA1724 
CHA1725 
CHA1726 
CHA1727 
xxxxxxxxCHA1728 
).  CHA1729 
xxxxxxxxCHA1730 
CHA1731 
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GO  TO  C1,2,3,4,5,6),NFLUX  CHA1732 

CONTINUE  CHA1733 

CHA1734 

NFLUX=1.    LINE  X=0  IS  TO  THE  LEFT  OF  LEFT  WAVE.  CHA1735 

CHA1736 

UX=UL  CHA1737 

PX=PL  CHA1738 

ROX=ROL  CHA1739 

ZX=ZL  CHA1740 

GX=GL  CHA1741 

DUDXIX=DUDXIL  CHA1742 

DPDXIX=DPDXIL  CHA1743 

DRDXIX  =  DRDXIL  CHA1744 

DZDXIX  =  DZDXIL  CHA1745 

DUDTX=-DPDXIL  CHA1746 

DR0DTX=-R0L**2*DUDXIL  CHA1747 

DPDTX=-GL**2*DUDXIL  CHA1748 

DRODTX=DRODTX-RAT*ROL*UL  CHA1749 

DPDTX=DR0DTX*CL**2  CHA1750 

DZDTX=0.  CHA1751 

GO  TO  9  CHA1752 

CONTINUE  CHA1753 

CHA1754 

NFLUX=6.    LINE  X=0  IS  TO  THE  RIGHT  OF  RIGHT  WAVE.                   CHA1755 

CHA1756 

UX=UR  CHA1757 

PX=PR  CHA1758 

ROX=ROR  CHA1759 

ZX=ZR  CHA1760 

GX=GR  CHA1761 

DUDXIX=DUDXIR  CHA1762 

DPDXIX=DPDXIR  CHA1763 

DRDXIX=DRDXIR  CHA1764 

DZDXIX=DZDXIR  CHA1765 

DUDTX=-DPDXIR  CHA1766 

DPDTX=-GR**2*DUDXIR  CHA1767 

DR0DTX=-R0R**2*DUDXIR  CHA1768 

DRODTX=DRODTX-RATxROR*UR  CHA1769 

DPDTX=DR0DTX*CRX*2  CHA1770 

DZDTX=0.  CHA1771 

GO  TO  9  CHA1772 

!    CONTINUE  CHA1773 

CHA1774 

NFLUX=2.    SONIC  CASE  (LEFT).  CHA1775 

CHA1776 

BETA0=(MU2*(UL/CL+G7) )**( 1 . D0/MU2)  CHA1777 

SQBO=DSQRT(BETAO)  CHA1778 

A1=DUDXIL+DPDXIL/GL  CHA1779 

A0=A1-(CL/(G15*SL) )*DSDXIL*( BETA0**G5-1 . DO )  CHA178  0 

EVER1=-((GAMA-1.D0)*UL+2.D0*CL)*(BETA0**G13-1.D0)/(GAMA-3.D0)  CHA1731 

EVER2=4.DO*CL*(BETAO**G14-l . DO )/( 3 . D0*GAMA-5 . DO)  CHA1782 

EVER=(EVER1+EVER2)*RAT/R0L  CHA1783 

A0=(A0+EVER)  CHA1784 

DPDAX=GL*BETAO*AO  CHA1785 

C0=MU2*(UL+G7*CL)  CHA1786 

IFCCO.LT.O.)  CALL  SOFC'FLUXE  2.   CO  NEGATIVE.')  CHA1787 

UX=C0  CHA1788 

ROX=GL*BETA0/C0  CHA1789 

ZX=ZL  CHA1790 

PX=ROX*C0**2/GAMA  CHA1791 

GX=ROX*C0  CHA1792 

DPDAX=DPDAX+RAT*UXXCO*SQBO  CHA17  93 

DUDBX=-CL*BETA0**(-1 .D0/G4)/G4  CHA1794 

DPDBX=PL*BETA0**MU2/G6  CHA1795 

DRODBX=ROL*BETA0**(-MU2)/G4  CHA17  96 

DSDAX=SQBO*DSDAL  CHA1797 

DZDAX=SQBO*DZDAL  CHA1798 

DRODAX=DPDAX/CO*X2-(ROX/(GAMA*SL))*DSDAX  CHA17  99 

DUDAX=A0  CHA1800 

DGDAX=0.5DO*GAMA*(PX*DRODAX+ROX*DPDAX)/GX  CHA1801 

GO  TO  9  CHA1802 

i           CONTINUE  CHA1803 
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FILE:  CHARGEP   FORTRAN   Al 


NFLUX=5 


SONIC  CASE  (RIGHT) 


BETA0=(MU2x(-UR/CR+G7))xx(l.D0/MU2) 

SQBO=DSQRT(BETAO) 

A1=DUDXIR-DPDXIR/GR 

A0=A1+(CR/(G15XSR))XDSDXIRX(BETA0X*G5-1.D0) 

EVERl=(-(GAMA-l.D0)xUR+2.D0XCR)X(BETA0XXG13-l.D0)/(GAMA-3 

EVER2=-4.D0XCRX(BETAOXXG14-l.DO)/(3.DOXGAMA-5.DO) 

EVER=(EVER1+EVER2)XRAT/R0R 

A0=(A0+EVER) 

DPDAX=-GRXBETA0xA0 

C0=MU2X(-UR+G7XCR) 

IFCCO.LT.O.)  CALL  SOFCFLUXE  5.   CO  NEGATIVE.') 

UX=-CO 

ROX=GR*BETAO/C0 

ZX  =  ZR 

PX=R0XXC0XX2/GAMA 

GX=ROXXCO 

DPDAX=DPDAX-RATXUXXCOXDSQRT(BETAO) 

DUDBX=CR*BETA0*X(-1.D0/G4)/G4 

DPDBX=PRXBETA0xxMU2/G6 

DRODBX=RORXBETA0XX(-MU2)/G4 

DSDAX=SQBOXDSDAR 

DZDAX=SQBOxDZDAR 

DRODAX=DPDAX/CO*X2-(ROX/(GAMA*SR))*DSDAX 
DUDAX=AO 

DGDAX=0.5D0XGAMAX(PXXDRODAX+ROXXDPDAX)/GX 

GO  TO  9 

CONTINUE 


DO) 


NFLUX=3 


LINE  X=0  IS  BETWEEN  THE  LEFT  WAVE  AND  THE  CONTACT. 


UX=USTAR 

PX=PSTAR 

ROX=RSTARL 

ZX  =  ZL 

GX=GSTARL 

DUDXIX=-DPIDT/GSTARLXX2 

DPDXIX=-DUIDT 

DUDXIX=DUDXIX-RATXUSTAR/RSTARL 

DZDXIX=DZDXIL 

DZDTX=0. 

IF  (  .NOT.HELEML)  GO  TO  32 

31  CONTINUE 
:   LEFT  SHOCK. 

DRDXIX=(RSTARL/WL)**2*(3.D0*DUIDT 

1  +DPIDT*(1.D0+3.D0X(WL/GSTARL)X*2)/WL 

2  +DUDXILXWLX((GL/WL)xx2+3.D0)+3.D0xDPDXIL 

3  +DRDXILX(WL/R0L)XX2) 
EVER1=ULXRSTARLXX2XRATX((GL/WL)XX2+1.DO)/(ROLXWL) 
EVER2=2 . DOXRSTARLXUSTARXRAT/WL 
DRDXIX=DRDXIX+EVER1+EVER2 
DR0DTX=-DUDXIXXR0XXX2 

GO  TO  33 

32  CONTINUE 
BETA=GSTARL/GL 
SQB=DSQRT(BETA) 
DPDA=ASTARLXGSTARL 

DPDA=GSTARLX(ASTARL+RATXUSTARXCSTARL/(GLX    SQB)) 
G41=l .D0/G4+0.5D0 

DR0DA=(DRDXIL-DPDXIL/(CLXCD)    *8ETAx*G41+DPDA/ (CSTARLXX2) 
DRDXIX=    DR0DA/SQB+DPIDT/CGSTARLXCSTARLXX2) 
DR0DA=DPDA/CSTARLX^2-(RSTARL/(GAMAXSL))XDSDASL 
DR0DTX=-DUDXIXXR0XXX2 
DRDXIX=DRODA/SQB+DRODTX/GSTARL 

33  CONTINUE 
DUDTX=DUIDT 
DPDTX=DPIDT 
GO  TO  9 

4     CONTINUE 


CHA1804 
CHA1805 
CHA1806 
CHA1807 
CHA1808 
CHA1809 
CHA1810 
CHA1811 
CHA1812 
CHA1813 
CHA1814 
CHA1815 
CHA1816 
CHA1817 
CHA1818 
CHA1819 
CHA1820 
CHA1821 
CHA1822 
CHA1823 
CHA1824 
CHA1825 
CHA1826 
CHA1827 
CHA1828 
CHA1829 
CHA1830 
CHA1831 
CHA1832 
CHA1833 
CHA1834 
CHA1835 
CHA1836 
CHA1837 
CHA1838 
CHA1839 
CHA1840 
CHA1841 
CHA1842 
CHA1843 
CHA1844 
CHA1845 
CHA1846 
CHA1847 
CHA1S48 
CHA1849 
CHA1850 
CHA1851 
CHA1852 
CHA1853 
CHA1354 
CHA1855 
CHA1856 
CHA1857 
CHA1353 
CHA1859 
CHA1860 
CHA1861 
CHA1862 
CHA1863 
CHA1864 
CHA1865 
CHA1866 
CHA1367 
CHA1368 
CHA1869 
CHA1870 
CHA1871 
CHA1872 
CHA1873 
CHA1874 
CHA1875 
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GEP   FORTRAN   Al 

NFLUX=4.    LINE  X=0  IS  BETWEEN  THE  CONTACT  AND  THE  RIGHT  WAVE.  CHA1876 

CHA1877 

DPDXIX=-DUIDT  CHA1878 

UX=USTAR  CHA1879 

PX=PSTAR  CHA1880 

ROX=RSTARR  CHA1881 

ZX=ZR  CHA1882 

GX=GSTARR  CHA1883 

DUDXIX=-DPIDT/GSTARRx*2  CHA1884 

DUDXIX=DUDXIX-RAT*USTAR/RSTARR  CHA1885 

DPDXIX=-DUIDT  CHA1886 

DZDXIX=DZDXIL  CHA1887 

DZDTX=0.  CHA1888 

IF  (.NOT.HELEMR)  GO  TO  42  CHA1889 

41  CONTINUE  CHA1890 
RIGHT  SHOCK  CHA1891 

DRDXIX=(RSTARR/WR)*X2*(3.*DUIDT  CHA1892 

1  -DPIDT*(1.D0+3.D0*(WR/GSTARR)**2)/WR  CHA1893 

2  -DUDXIRXWRX((GR/WR)X*2+3.D0)  +  3.D0*DPDXIR  CHA1894 

3  +DRDXIRX(WR/R0R)XX2)  CHA1895 

EVER1=UR*RSTARR**2*RAT*((GR/WR)**2+1 . D0)/(RORXWR)  CHA1896 

EVER2=2.D0xRSTARRxUSTARXRAT/WR  CHA1897 

DRDXIX=DRDXIX-EVER1-EVER2  CHA1898 

DR0DTX=-DUDXIXXR0XXX2  CHA1899 

GO  TO  43  CHA1900 

42  CONTINUE  CHA1901 
RIGHT  RAREFACTION  CHA1902 

BETA=GSTARR/GR  CHA1903 

SQB=DSQRT(BETA)  CHA1904 

DPDA=-ASTARRXGSTARR  CHA1905 

DPDA=-GSTARRX(ASTARR+RATXUSTARXCSTARR/(GRX    SQB))  CHA1906 

G41=1.D0/G4+0.5D0  CHA1907 

DRODA=(DRDXIR-DPDXIR/(CR*CR))    XBETAX*G41+DPDA/(CSTARRXX2)  CHA19  08 

DRDXIX=    DR0DA/SQB-DPIDT/(GSTARR*CSTARR**2)  CHA1909 

DR0DA=DPDA/CSTARRX*2-(RSTARR/(GAMAXSR))XDSDASR  CHA1910 

DR0DTX=-DUDXIXXR0XX*2  CHA1911 

DRDXIX=DRODA/SQB-DRODTX/GSTARR  CHA1912 

43  CONTINUE  CHA1913 
DUDTX=DUIDT  CHA1914 
DPDTX=DPIDT  CHA1915 
GO  TO  9  CHA1916 

9  CONTINUE  CHA1917 

XXXXX^XXX^XXX^^)f^^^X^^¥^XX^X5«5«X5(^X^5(5(XXXX5(X^5(^^5(^X^3(5(3(X5(X^X3«XXX^X^5(^5€^XCHA1918 

FLUXES  CENTERED  AT  TIME  TCN+1/2)  AT  EULERIAN  POINT  X=X(I).  CHA1919 

^XXX^^)(*^^^X^^^^^^XX^^X^X^^X5(XX^XX^^XX3«XXX^^^^XXXXX3(5(XXXXXXXXX^^XX^XXX)(CHA1920 

FIl=ROXXUX  CHA1921 

FI2=R0X*UX**2+PX  CHA1922 

FI2=FI2-PX  CHA1923 

FI3=UX*(G12*PX+0.5DO*ROXXUXX*2)  CHA1924 

FI4=ZXXR0XXUX  CHA1925 

FI3=FI3+QDETXFI4  CHA1926 

ROUOO=ROX*UX  CHA1927 

GO  T0(10,20,30,40,50,60),  NFLUX  CHA1923 

10  CONTINUE  CHA1929 
60    CONTINUE  CHA1930 

DFDXI1=DRDXIX*UX+R0X*DUDXIX  CHA1931 

DFDXI2=DRDXIXXUX**2+2.DO*ROX*UXXDUDXIX+DPDXIX  CHA1932 

DFDXI2=DFDXI2-DPDXIX  CHA1933 

DFDXI3=DUDXIXX(G12*PX+0.5D0xROXXUXX*2)  CHA19  34 

1      +UXx(G12xDPDXIX+0.5D0xDRDXIXXUXX*2+ROXxUXxDUDXIX)  CHA1935 

DFDXI4=ZXXDFDXI1+R0XXUXXDZDXIX  CHA19  36 

DFDXI3=DFDXI3+QDET*DFDXI4  CHA19  37 

DFIDTl=DRODTXXUX+ROXXDUDTX  CHA19  38 

DFIDT2=DR0DTXXUX**2+2.D0XR0XXUX*DUDTX+DPDTX  CHA19  39 

DFIDT2=DFIDT2-DPDTX  CHA1940 

DFIDT3=DUDTXX(G12XPX+0.5DOXROXXUXX^2)  CHA1941 

1      +UXX(G12XDPDTX+0.5D0XDR0DTXXUXX^2+R0XXUXXDUDTX)  CHA1942 

DFIDT4=ZXXDFIDT1+R0XXZXXDZDTX  CHA1943 

DFIDT3=DFIDT3+QDET*DFIDT4  CHA1944 

FIDOT1=-ROUOOXDFDXI1+DFIDT1  CHA1945 

FIDOT2=-ROUOOXDFDXI2+DFIDT2  CHA1946 

FIDOT3=-ROU00XDFDXI3+DFIDT3  CHA1947 
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FILEs  CHARGEP   FORTRAN   Al 


FIDOT4=-ROUOO*DFDXI4+DFIDT4  CHA1948 

UXDOT=-ROUOO*DUDXIX+DUDTX  CHA1949 

PXDOT=-ROUOO*DPDXIX+DPDTX  CHA1950 

ROXDOT=-ROUOO*DRDXIX+DRODTX  CHA1951 

ZXDOT=-ROUOO*DZDXIX+DZDTX  CHA1952 

FIH1=FI1+DT2*FID0T1  CHA1953 

FIH2=FI2+DT2*FID0T2  CHA1954 

GIH=PX+DT2*PXD0T  CHA1955 

FIH3=FI3+DT2*FID0T3  CHA1956 

FIH4=FI4+DT2*FID0T4  CHA1957 

UXN=UX+DT*UXDOT  CHA1958 

PXN=PX+DT*PXDOT  CHA1959 

ROXN=ROX+DT*ROXDOT  CHA1960 

ZXN=ZX+DT*ZXDOT  CHA1961 

IFCZXN.LT.O.)  ZXN=0.  CHA1962 

GO  TO  90  CHA1963 

20    CONTINUE  CHA1964 

EVO=GL*DSQRT(BETAO)  CHA1965 

201  CONTINUE  CHA1966 

DFIDA1=DR0DAX*UX+R0X*DUDAX  CHA1967 

DFIDA2=DR0DAX*UX**2+2 .  D0*ROX*UX*DUDAX+DPDAX  CHA1968 

DFIDA2=DFIDA2-DPDAX  CHA1969 

DFIDA3=DUDAX*(G12*PX+0.5D0*ROX*UX**2)  CHA1970 

1      +UX*(G12*DPDAX+0.5D0*DRODAX*UX**2+ROX*UX*DUDAX)  CHA1971 

DFIDA4=ZX*DFIDA1+R0X*UX*DZDAX  CHA197  2 

FIDOT1=-EVO*DFIDA1  CHA1973 

FIDOT2=-EV0*DFIDA2  CHA1974 

FIDOT3=-EV0*DFIDA3  CHA1975 

FIDOT4=-EV0*DFIDA4  CHA1976 

FIH1=FI1+DT2*FID0T1  CHA1977 

FIH2=FI2+DT2*FID0T2  CHA1978 

FIH3=FI3+DT2*FID0T3  CHA1979 

FIH4=FI4+DT2*FID0T4  CHA1980 

GA=DGDAX  CHA1981 

IF(NFLUX.EQ.5)GA=-GA  CHA1982 

DROUA=UX*DRODAX+ROX*DUDAX  CHA1983 

BETAPR=0.5D0*DSQRT(BETA0)*(GA-DROUA)  CHA1984 

FIH2=FIH2-DPDBX*BETAPR*DT2  CHA1985 

UXDOT=-EV0*DUDAX+BETAPR*DUDBX  CHA1986 

PXDOT=-EV0*DPDAX+BETAPR*DPDBX  CHA1987 

GIH=PX+DT2*PXD0T  CHA1988 

ROXDOT=-EV0*DRODAX+BETAPR*DRODBX  CHA1989 

ZXDOT=-EVO*DZDAX  CHA1990 

UXN=UX+DT*UXDOT  CHA1991 

PXN=PX+DT*PXDOT  CHA1992 

ROXN=ROX+DT*ROXDOT  CHA1993 

ZXN=ZX+DT*ROXDOT  CHA1994 

IFCZXN.LT.O.)  ZXN=0.  CHA1995 

GO  TO  90  CHA1996 

50    CONTINUE  CHA1997 

EVO=-GR*DSQRT(BETAO)  CHA1998 

GO  TO  201  CHA1999 

30    CONTINUE  CHA2000 

40    CONTINUE  CHA2001 

GO  TO  60  CHA2002 

90    CONTINUE  CHA2003 

RETURN  CHA2004 

END  CHA2005 
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Figure  A-l. 


Piecewise    Linear   Distribution   of  Flow   Variables   in   Cells 


53 


gAGEP/icr/o/v/   BRAkjcH 


SHOCK    BP/^AICH 


Vl  U* 


u* 


Figure  A-2.  Intersection   of  Right   and   Left  Adiabats   for  Solving   Riemann   Problem 
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Figure  A-3. 


Wave   Diagram   Representing   Solution   to    Riemann   Problem 
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APPENDIX  B.      Code  for  Re-Normalizing  the  Air  Impulse 


7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 


.9,1., 


221, .178, .149, 
0136, .0113, .0095/ 


IMPLICIT  REAL*8(A-H,0-Z) 
C   CODE  RENORM  --  C   TRANSFORMATION  OF  TOTAL  REFLECTED  IMPULSE  FROM 
C   BAKER'S  CHART  TO  SPACE-NORMALIZED  VALUES. 
C   DATA  FROM  FIG.  6.3  (SUPPLEMENT)  IN  BAKER'S  BOOK  "EXPLOSIONS  IN  AIR" 

REAL*4  RB,IB,RS,IS,ISBARE 

DIMENSION  RB(21),IB(21) 

DIMENSION  RS(21),IS(21),ISBARE(21) 

DATA  RB/. 05, .06, .07, .08, .09, .1, .2, .3, .4, .5, .6, .7, .8, 
1        2. ,3. ,4. ,5. ,6. ,7./ 

DATA  IB/4.4,3.06,2.30,1.83,1.50,1.27, .457, .293, 
1         .128, .113, .099, .0885, .0376, .0236, .0173, 

PAI=4.D0*DATAN(1.D0) 

G=1.4D0 

PA=0.1D0 

RH0A=1.3D0 

RHO0=1800.D0 

Q0=4.D0 

BETA=DSQRT(RHOA/RHO0)*(PA/(RHO0*Q0))**(l.D0/6.D0) 

GOREM=(3.D0/DSQRT(2.D0*G))*(4.D0*PAI/3.D0)*x(l.D0/3.D0) 

BETA=BETA*GOREM 

DELTA=(  (4.D0*PAI/3.D0)*(RHO0XQ0/PA)  )**C1 . D0/3.D0) 
PRINT  11,  BETA, DELTA 
11    FORMATC/1X, 'RESULTS  WITH   BETA, DELTA= • , 2D16 . 7// 


IX,  » 
i 


N 


RS 


RB 


IS 


IB 
',2X,' 


,2X, 

ISBARE 


'/) 


19 
20 
21 
22 

23 
24 
25 

26 


DO  1  N=l,21 

RS(N)=RB(N)XDELTA 

IS(N)=IB(N)*BETA 

ISBARE(N)=1.D0/RS(N)XX2 

PRINT  2,  N,RB(N),IB(N),RS(N),IS(N),ISBARE(N) 

F0RMAT(1X,I4,2E12.4,2X,2E12.4,2X,E12.4) 

CONTINUE 

END 


RESULTS  WITH   BETA,DELTA= 
N      RB  IB 


0.1204163D-01 
RS 


0.6706157D+02 

IS  ISBARE 


1 
2 
3 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 


5000E-01 
6000E-01 
7000E-01 
8000E-01 
9000E-01 
1000E+00 
2000E+00 
3000E+00 
4000E+00 
5000E+00 
6000E+00 
0.7000E+00 
0.3000E+00 
9000E+00 
1000E+01 
2000E+01 
3000E+01 
4000E+01 
5000E+01 
6000E+01 
7000E+01 


0.4400E+01 
0.3060E+01 
0.2300E+01 
0.1830E+01 
0.1500E+01 
0.1270E+01 
0.4570E+00 
0.2930E+00 
2210E+00 
1780E+00 
1490E+00 
1280E+00 
1130E+00 
9900E-01 
0.8850E-01 
0.3760E-01 
0.2360E-01 
0.1730E-01 
0.1360E-01 
0.1130E-01 
0.9500E-02 


3353E+01 
4024E+01 
4694E+01 
5365E+01 
6036E+01 
6706E+01 
1341E+02 
2012E+02 
2682E+02 
3353E+02 
4024E+02 
4694E+02 
5365E+02 
6036E+02 
6706E+02 
1341E+03 
2012E+03 
2682E+03 
3353E+03 
4024E+03 


0.4694E+03 


0.5298E-01 
0.3685E-01 
0.2770E-01 
0.2204E-01 
0.1806E-01 
1529E-01 
5503E-02 
3528E-02 
2661E-02 
0.2143E-02 
0.1794E-02 
0.1541E-02 
0.1361E-02 
0.1192E-02 
0.1066E-02 
0.4528E-03 
0.2842E-03 
0.2083E-03 
0.1638E-03 
0.1361E-03 
0.1144E-03 


8894E-01 
6177E-01 
4533E-01 
3474E-01 
2745E-01 
2224E-01 
5559E-02 
2471E-02 
1390E-02 
8894E-03 
6177E-03 
4538E-03 
0.3474E-03 
0.2745E-03 
2224E-03 
5559E-04 
2471E-04 
1390E-04 
8894E-05 
0.6177E-05 
0.4538E-05 
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